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ABSTRACT 


Procedures  for  designing  feedback  controllers  for 
maintaining  a  prescribed  level  of  dissolved  oxygen  in  a 
single  reach  of  a  polluted  river  system  are  described  in 
this  thesis.  A  single  control  scheme  in  which  artificial 
instream  aeration  is  used  as  the  control  and  a  two-control 
scheme  using  aeration  and  effluent  discharge  as  the  two 
controls  are  developed.  A  discrete-time  model  of  the 
polluted  river  system  used  in  the  controller  design  is 
derived  first  starting  from  an  available  continuous  time 
model . 

Unlike  continuous-time  controllers  developed  by 
others  previously,  the  advantage  of  the  discrete-time 
controllers  is  that  the  output  DO  (dissolved  oxygen)  does 
not  have  to  be  monitored  at  all  times.  This  is  very 
attractive  from  a  practical  point  of  view. 

The  effect  of  seasonal  variations  in  the 
temperature  of  the  water  on  the  performance  of  the  digital 
controllers  developed  in  this  thesis  as  well  as  the 
continuous  controllers  reported  in  the  literature  is  also 
investigated  by  means  of  computer  simulations. 
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CHAPTER  1 


INTRODUCTION 

1 . 1  The  Environment 

Man  has  not  paid  much  attention  to  his  environ¬ 
ment  until  recently.  In  as  much  as  growth,  achieved  by 
increased  industrialization  and  urbanization,  is  desired 
in  our  societies,  it  is  essential  that  steps  are  taken 
to  alleviate  the  problems  that  accompany  such  growth.  One 
of  these  problems  is  environmental  pollution.  The  effect 
of  man's  activities  on  environment  in  or  near  industrial 
and  urban  centers  is  just  being  known.  The  increase  in 
the  rate  of  discharge  of  industrial  and  urban  waste  into 
the  environment  has  not  been  matched  by  an  equivalent 
increase  in  the  assimilative  capacity  of  the  environment, 
with  the  consequent  degradation  of  the  basic  requirements 
for  life  itself;  (viz)  air  and  water. 

1 . 2  Sources  of  Environmental  Pollution 

It  is  evident  that  the  soil  which  is  the  main 
source  of  man's  food  has  been  turned  into  a  dumping  ground 
for  pollutants  such  as  solid  refuse  from  both  urban  and 
rural  sources,  chemical  agents  used  in  the  production  of 
food  and  other  industrial  processes  and  particles  settling 
out  of  the  air  as  a  result  of  air  pollution.  It  is  now 
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known  that  some  of  these  pollutants  cause  harmful  residues 
to  be  concentrated  in  the  food  chain  [ 9 ] 1 . 

The  pollution  of  the  air  is  caused  by  the  burning 
of  fuels  for  transportation,  heat  and  power  and  the  emission 
of  sulphur  dioxide,  particle  matter,  carbon  monoxide, 
photo-chemical  oxidants,  hydrocarbons  and  nitrogen  dioxide 
into  the  air  by  industrial  processes  that  give  off  these 
gases  as  waste  products.  Environmentalists  speculate  that 
the  oxygen-carbon  dioxide  balance  might  be  upset  due  to 
these  pollutants  [9].  In  the  case  of  water  pollution,  the 
use  of  lakes,  rivers  and  the  oceans  by  municipalities  and 
industries  for  the  dilution  of  their  liquid  wastes  and  as 
dumping  ground  for  their  solid  wastes  has  made  these  places 
unsafe  for  recreational  activities  such  as  swimming  and 
most  important  as  a  source  of  food;  fish. 

1 . 3  Water  Pollution 

There  are  many  factors  to  be  considered  in  assessing 
the  quality  of  water  in  a  river  or  stream.  These  are: 

Physical  factor:  The  color  of  the  water. 

Chemical  factors:  Amount  of  the  various  chemical  compounds 
present,  e.g.  nitrates,  chlorides,  phosphates,  etc. 

Biological  factors:  The  composition  of  life  in  the  river. 

It  is  a  fact  that  in  a  normally  healthy  river,  there  is  a 

1  The  number  in  square  brackets  refer  to  the  references  listed 
under  bibliography  at  the  end  of  this  thesis. 
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biodynamic  cycle  which  results  in  a  balance  between  the 
plant  and  animal  life.  Ruth  Patrick  [13]  has  developed 
a  system  for  assessing  the  degree  of  pollution  of  a  river 
from  an  analysis  of  the  various  groups  of  organisms 
present  in  the  river. 

Biochemical  factors:  The  biochemical  oxygen  demand,  BOD 
and  the  dissolved  oxygen,  DO.  The  biochemical  oxygen 
demand  is  defined  as  the  amount  of  oxygen  required  by 
bacteria  to  decompose  a  given  amount  of  organic  waste  under 
aerobic*  conditions  at  a  given  temperature.  Dissolved 
oxygen,  DO,  is  defined  as  the  amount  of  oxygen  absorbed  by 
a  river  at  a  given  temperature. 

Of  all  the  factors  stated  above  BOD  and  DO  are 
generally  accepted  as  the  primary  indicators  of  water 
quality.  The  reasons  for  this  will  be  clear  from  a 
discussion  of  the  pollution  of  a  river  by  organic  wastes 
as  outlined  by  Brinley  [13] . 

Organic  waste  dumped  into  the  river  is  actually 
broken  down  by  bacteria  present  in  the  river.  This 
decomposition  is  characterized  by  three  phases.  These  are: 
(i)  Active  phase 

In  this  phase  there  is  increased  bacterial  activity 
resulting  in  increased  respiration  and  therefore 
increased  oxygen  consumption.  The  biochemical  oxygen 

■k 

Reaction  occurring  only  in  the  presence  of  oxygen. 
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demand  is  high  and  thus  a  low  level  of  dissolved 
oxygen  results.  The  biodynamic  cycle  is  disrupted 
and  organisms  that  require  a  lot  of  oxygen  in  order 
to  survive  become  extinct. 

(ii)  Intermediate  phase 

Since  most  of  the  waste  has  already  been  decomposed 
in  the  active  phase,  the  amount  of  energy  required 
by  the  bacteria  to  break  down  the  remaining  waste  is 
greatly  reduced  in  this  phase.  Dissolved  oxygen  level 
rises  as  a  result  of  this  and  some  living  organisms 
start  to  return  to  the  river. 

(iii)  Recovery  phase 

As  the  waste  is  progressively  decomposed,  the  water 
in  the  river  becomes  clearer,  and  the  biodynamic  cycle 
comes  into  effect  once  again. 

Thus  it  can  be  seen  that  low  levels  of  biochemical 
oxygen  demand  and  high  levels  of  dissolved  oxygen  are 
characteristics  of  a  river  free  from  organic  waste  pollution. 

1 . 4  Application  of  Automatic  Control  Theory  to  Water 

Pollution 

It  is  evident  from  the  last  section  that  it  takes 
some  time  for  the  river  to  recover  from  the  effects  of 
organic  waste  pollution.  Thus  the  greater  the  rate  of 
waste  discharged  into  the  river,  the  less  the  river  system 
is  able  to  assimilate  it  .  In  the  past  few  years  some 
researchers  have  focused  their  attention  on  the  engineering 
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aspects  of  improving  and  controlling  the  water  quality  in 
rivers  which  are  subjected  to  organic  waste  dumping.  Some 
of  them  [8,  22]  have  formulated  the  water  pollution  problem 
as  an  optimization  problem.  Not  only  is  such  a  formulation 
complicated  but  in  most  cases  the  solutions  arrived  at 
are  not  easy  to  implement.  Young  and  Beck  [24]  have  used 
a  simpler  approach  to  the  problem.  They  have  used  the  pole 
assignment  technique  to  design  a  controller  that  will 
maintain  a  prescribed  level  of  DO  in  a  single  reach  of  a 
river  system  using  complete  state  feedback.  Ramar  and 
Gourishankar  [18,19]  have  gone  further  and  designed 
controllers  using  output  feedback  (DO)  only  for  pole  assign¬ 
ment  . 

They  have  also  used  the  minimization  of  a  sensi¬ 
tivity  function  to  obtain  a  better  response  than  those 
obtained  by  Young  and  Beck.  Kudva  and  Gourishankar  [11] 
in  their  work  in  this  field  have  designed  an  observer  to 
estimate  the  output  BOD  whose  exact  and  continuous  measure¬ 
ment  is  difficult  and  have  been  able  to  show  that  an  increased 
effluent  discharge  on  the  whole  over  a  period  of  time  can 
be  achieved  by  their  method. 

In  all  these  cases  the  section  of  the  river  whose 
water  quality  is  to  be  maintained  has  been  modelled  as  a 
continuous  system.  This  means  that  the  output  DO  has  to 
be  monitored  or  the  BOD  estimated  all  the  time.  A  more 
practical  approach  would  be  to  monitor  the  output  DO  only 
at  certain  predetermined  instants  of  time  and  use  this 
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information  to  control  the  output  DO  during  the  period 
between  samplings.  Such  a  discrete  time  controller  would 
be  easier  to  implement,  particularly  in  this  age  of 
digital  computers. 

1 • 5  Scope  of  this  Thesis 

The  work  reported  in  this  thesis  is  primarily 
concerned  with  the  design  of  a  direct  digital  controller 
(DDC)  to  maintain  a  prescribed  level  of  DO  in  a  single 

reach  of  a  river  system  subjected  to  dumping  of  organic 
waste . 

A  comparative  study  is  also  made  between  the 
performance  of  the  digital  controller  and  that  of  the 
continuous  time  controllers  designed  by  Ramar  and  Gourishankar , 
or  Kudva  and  Gourishankar  [18,  19,  10]. 

In  Chapter  Two  of  this  thesis,  a  discrete-time 
state  space  model  for  a  single  reach  of  a  river  system  is 
derived.  Two  cases  are  treated.  One  with  artificial  aeration 
as  the  only  control  and  the  other  with  aeration  plus  effluent 
discharge  as  the  two  controls. 

Chapter  Three  describes  the  design  of  the  digital 
controllers  for  the  two  cases  mentioned  above,  using  the 
pole-assignment  and  output  (DO)  feedback. 

The  effect  of  seasonal  temperature  variation  on 
the  performance  of  the  continuous  and  digital  controllers  in 
maintaining  water  quality  is  investigated  in  Chapter  Four. 
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comparisons  are  made 
for  further  research 


Five,  the  results  are  summarized, 
and  conclusions  drawn.  Suggestions 
are  also  given. 


CHAPTER  2 


DISCRETE-TIME  STATE  SPACE  MODEL 


OF  A  RIVER  SYSTEM 


2.1  Continuous-Time  Model 


As  mentioned  in  the  previous  chapter  the  most 
commonly  accepted  parameters  for  assessing  water  quality 
are  the  concentration  of  dissolved  oxygen  (DO)  and  the 
biochemical  oxygen  demand  (BOD)  in  a  section  of  the  river. 

A  number  of  different  models  have  been  derived  over  the 
past  few  years  starting  from  the  general  diffusion  equation 
and  using  BOD  and  DO  as  variables  [1,  4,  6,  21,  22]. 

One  such  model  used  by  many  researchers  is  given 
in  equation  (2-1)  below. 


DO 
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dt 
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(2-1) 


Where 


x^(t)  is  the  output  (downstream)  DO  in  mg/& 
x2(t)  is  the  output  BOD  in  mg/£ . 

L^(t)  is  the  input  (upstream)  BOD  in  mg/£ . 
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C^(t)  is  the  input  DO  in  mg/£. 

C  is  the  saturation  value  of  DO  in  mg/£. 
s 

L„  is  the  BOD  in  the  effluent  in  mg/£. 

E 

is  the  DO  in  the  effluent  in  mg/£. 

E 

is  the  mean  rate  of  addition  of  BOD  by  local 
A 

runoff . 

D  is  the  net  rate  of  removal  of  DO  by  the  combined 

action  of  photosynthesis  and  respiration. 

V  is  the  mean  volume  of  the  water  held  in  the 
m 

reach  in  m3 . 

Q  is  the  volumetric  flow  rate  in  m3/day. 

Q„  is  the  volumetric  flow  rate  of  the  effluent 
E 

in  m3/day. 

a^  is  the  reaeration  rate  in  day-1. 
a2  is  the  BOD  decay  rate  in  day-1. 

Equation  (2-1)  can  be  rewritten  in  vector 
matrix  form  as 
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where  x1 (t)  and  x2 (t)  represent  the  differentiation  of 
x1(t)  and  x 2 ( t )  with  respect  to  time. 
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The  continuous-time  model  given  in  equation  (2-2) 
does  not  include  the  control  term.  Introducing  this  term, 
the  continuous-time  model  with  artificial  aeration  as  the 
only  control  becomes 


Xx  (t) 

Q  +  Q 


+ 


E. 


V 


m 


0 


-  a, 


Q  +  Q, 


(a2  +  “ V 


m  — 1 


x1  (t) 


x2(t) 


+ 


Q 

V 

m 

0 

C±(t) 

+ 

qecf 

anC  -  Dd  +  fy  E 
Is  B  V 

m 

0 

lO 

l 

Li(t) 

qflf 

L  +  E  E 

A  V 

m 

+ 


1 

0 


UA(t> 


(2-3) 


Where  UA(t)  is  the  aeration  control. 

Another  method  of  controlling  the  output  DO  in 
the  river  is  to  use  not  only  artificial  aeration  but  also 
controlled  dumping  of  effluents  as  the  second  control.  By 
effluent  discharge  control  is  meant  the  regulation  of  the 
amount  of  waste  discharged  into  the  river.  A  good  control 
parameter  would  be  the  deviation  of  the  quantity  of  waste 
discharged  Q_  from  a  nominal  value  Q  ,  i.e. 

hi  hi 


qe  qe  +  6qe 


(2-4) 
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Substituting  equation  (2-4)  into  the  continuous¬ 
time  model  represented  by  equation  (2-1)  and  adding  the 
aeration  control  term  we  have 
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Let  U„  = 
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(2-5a) 


m 


Substituting  equation  (2-5a)  into  (2-5) ,  equation 


(2-5)  can  be  written  in  vector  matrix  form  as 


X1 

x2  (t) 


r  Q  +  qp 

— i 

P  — I 

-(al+  V  +UE> 

a2 

X1  (t) 

m 

q  +  qe 

0 

(a2  +  v 

+  UE> 

x2(t) 

m 

_  _ 

Q_C 

0-  0 
Vm 

Ci(t) 

+ 

alCs-DB+  V 

m 

Q  L_ 

0  — 

V 

Li(t) 

L  +  E  E 

A  V 

m 

m 

1 

CE 

uA(t) 

0 

le 

uE(t) 

(2-6) 
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where  U^(t)  is  the  aeration  control  and 

Ug ( t)  is  the  effluent  discharge  control. 

Equation  (2-6)  is  the  continuous-time  model  for  two  controls. 

2 . 2  Discrete-Time  Models 

The  river  system  to  be  controlled  is  inherently 
a  continuous  process.  In  this  case  the  digital  signal 
obtained  by  sampling  the  error  between  the  actual  output 
and  the  desired  output  will  first  be  processed  by  the 
digital  controller;  the  output  of  the  digital  controller 
f*(t) ,  which  is  a  pulse  train  is  then  smoothed  by  a  data- 
reconstruction  device  before  the  signal  is  applied  to  the 
continuous  process,  the  river  system. 

Figure  2.1  shows  the  schematic  diagram  of  this 
control  scheme.  The  data-reconstruction  device  is  a  zero- 
order  hold. 


Other  Inputs 


Figure  2.1  Closed  Loop  Discrete-Time  System 


In  continuous-time  systems,  the  state  equations  are 


a  set  of  first-order  differential  equations.  For  discrete- 
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time  systems  however,  the  state  equations  are  a  set  of 
first-order  difference  equations.  A  set  of  first-order 
difference  equations  can  be  obtained  from  the  set  of  first- 
order  differential  equations  by  taking  a  look  at  the 
combination  of  the  zero-order  hold  devices  and  the  contin¬ 


uous  system  shown  in  the  dotted  box  in  figure  2.1. 

Let  the  other  inputs  be  C^(t),  L^(t)  and  0. 

h  (t)  =  f  (KT)  =  U  (KT)  for  KT  £  t  £  (K  +  1)T  (2-7) 

A 

Since  h(t)  is  the  input  to  the  continuous-time  model. 

Also  m1(t)  =  C^KT)  (2-8) 

m2(t)  =  Li(KT)  (2-9) 

m3(t)  =  0 (KT)  (2-10) 


for  KT  <  t  <  (K  +  1)T 


where 


0  = 


an  C 
1  s 


Q  +  Q 


Q  +  Q, 


la  + 


V 


m 


2 -2a  Artificial  aeration  control  only 

Taking  Laplace  transforms  of  equation  (2-3)  we 


have 


s  X  (s)  -  x(tQ) 


Q  +  Qt 


(a1  + 


V. 


-  a, 


m 


Q  +  Q, 


(a2  +  V 


m 


X(s)  + 


0 
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—  — 

Q  0 

V 

m 

Ci(s) 

+  i  + 

1 

0  4 

Li  (s) 

s 

0 

_  _ 

_ 

___  ___ 

where 


xl(tJ 

x1  (s) 

i(tO>  = 

X(s)  = 

_X2(tol 

x2  (s) 

and 


0  = 


Q  +  Q 


aiCS  "  DB  + 


E 


V 


m 


Q  +  Q, 


la  + 


V 


m 


Let  A 


Q  +  Qt 


(a1  + 


V 


-  a 


m 


2 

Q  +  Q, 


0 


(a2  + 


v_ 


m 


Solving  for  X(s)  we 
X (s)  =  (si  - 


(Si  -  A)-1  = 


have 


x(t  )  + 

rQ 

V 

m 

0 

Ci(s) 

+  i 

0 

Q 

V 

m 

L±  (s) 

b 

1 

+ 

_0_ 

UA(s) 

1  a2 
(s  +  Kx)  (s  +  K1)  (s  +  K2) 

1 

s  +  K2 


s) 

(2-11) 


(2-12) 


0 


(2-13) 
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Q  +  Q, 


where  =  a^  +  — y~ 


m 


Q  +  Q, 


K. 


a2  + 


V 


m 


Also  taking  Laplace  transforms  of  equations  (2-7),  (2-8), 

(2-9)  and  (2-10)  we  have 


H  ( s )  = 


ua(kt) 


C. (KT) 

M  (s)  =  -i- - 

1  s 


L. (KT) 
M  (S)  =  - 

z  s 


M3(s)  = 


0  (KT) 


(2-14) 


Substituting  equations  (2-13)  and  (2-14)  into  (2-12)  and 
taking  inverse  Laplace  transforms,  the  resulting  equation 
can  be  written  in  vector  matrix  form  as 


x ( t )  =  $(t  -  tQ)  x(tQ)  +  B ( t  -  tQ)  w (KT)  +  D ( t  -  tQ) 


+  E(t  -  tQ)UA(KT) 


(2-15) 


where 


*(t  -  to* 


-K^t  -  tD) 


-  a. 


a2  al 


— -K1 (t 
e 


to) 


0 


e'V^V 


e~K2  ^t_to) 


(2-16) 
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B(t-tQ) 


2 'v  '  K,e 
m— [1 — 2 


-K]_  (t-tQ) 


K1K2 


a2  al 


+ 


Kie-K2(t-to)] 

a2“al 

(|_)  [i_e-K2(t-to)] 
vm 

K~ 


(2-17) 


D ( t-t  ) 

O 


0. 


0 


Kxe 


-K2(t-tQ) 


-] 


a2  al 


0X2  [1-e  K2 (t  to) ] 
_ 


-Kl(t-to) 


a2  al 


+ 


(2-18) 


E(t-to)  = 


l-e’Kl(t"to) 


K- 


0 


(2-19) 


w  (KT) 


C. (KT) 

l 


L± (KT) 


(2-20) 


0 


xi 


qfcf 

=  alCS  -  °B  +  — 


0 


Q  +  Q, 


m 


x2 


=  LA  + 


Q  +  Q 


K. 


al  + 


E 


V 
m 

Q  +  Q, 


V 


K, 


a2  + 


V 


m 


m 


Let  t  =  KT,  t  =  (K+l )  T ,  therefore  t-t  =  T  which 
o  ° 

is  the  sampling  interval.  Substituting  for  t,  tQ  and  t  -  tQ 


b 
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in  equation  (2-15) ,  we  have 


x  (K  +  1 )  T  =  $(T)  x  (KT)  +  B  (T)  w  (KT)  +  D  (T) 

+  E (T)  Ua(KT)  (2-21) 

where  $  (T)  ,  B(T),  D(T)  AND  E  (T)  are  given  by  equations 
(2-16) ,  (2-17) ,  (2-18)  and  (2-19)  with  t  -  t  replaced  by  T. 

Equation  (2-21)  is  a  set  of  first-order  difference 
equations  and  therefore  is  the  discrete  time  model  of  the 
river  system  with  aeration  control  only. 

2 -2b  Artificial  aeration  plus  effluent  discharge  controls 
Equation  (2-6)  can  be  written  as 


x(t)  =  A(t)  x(t)  +  bw(t)  +  d  +  eu(t) 


(2-22) 


where 


A  ( t)  = 


Q  +  Q 


(ai +  — 


E 


m 


+  uE(t)) 


-  a. 


0 


Q  +  Q, 


^a2  +  V 


m 


+  UE(t) 


b  = 


_Q_ 

V 


m 


0 


0 

Q_ 

V 

m 


d  = 


QfCf 

alCs  "  °B  +  ~ 

m 


^ELE 

la  +  — 

m 
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ce" 

°A(t) 

C±(t) 

U(t)  = 

w(t) 

_°  V 

_UE(tl 

Li(t) 

Since  A(t)  is  a  time  varying  matrix,  the  derivation 

of  the  state  transition  matrix,  using  the  Laplace  transforms 

technique  becomes  complicated  and  in  fact  impossible  since 

the  control  function  U_(t)  is  still  to  be  found.  In  order 

E 

to  simplify  the  derivation  of  the  discrete-time  model  for 
the  two  controls  case,  using  the  Laplace  transform  technique, 
it  will  be  assumed  that,  since  U^t)  varies  very  slowly, 

U„(t)  is  constant  over  a  sample  interval  KT  <  t  <  (K  +  1)T. 
Therefore  equation  (2-22)  can  be  written  as 


x(t)  =  A  (KT)  x(t)  +bw(t)  +  d  +  eu(t) 


(2-23) 


for  KT  <  t  <  (K  +  1)T 


where 


A  (KT)  =  A  ( t ) 


and  K  =  0 ,  1 ,  2 


t  =  KT 


A (KT)  is  now  a  constant  matrix  between  sampling 
instants,  hence  a  set  of  difference  equations  can  be  obtained 
using  the  Laplace  transform  technique  as  in  section  2-2a. 

The  discrete-time  model  is 


x (K  +  1 ) T  =  $(K,T)  X (KT)  +  B(K,T)  w(KT)  +  D(K,T) 

+  E  (K , T)  u  (KT)  K  =  0,  1,  2  -  (2-24) 
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whei*e 


<J>  (K,T)  = 


rciT 


a2 


a2  al 


[e 


~C1T  _  “c2t- 


B  (K , T)  = 


0 


{£-)  (1  ~  e  ClT) 
m 


rc2T 


-a. 


C1C2 


[1  - 


-C-|T  -CnT 

c~e  x  c,  e  ^ 

— - +  — - ] 

a2  al  a2  al 


Q 


-c2T, 


0 


(f~)  (1  -  e 
vm 


D  (K ,  T)  = 


4>x  (1  -  e  clT 

X1 


)  i  -CiT  -CiT 

a^d>x0  c^e  -1  c^e  z 


2  ■  2 
C1C2 


[1  - 


a2  al  a2  al 


] 


x 


-  [1  -  e  c2T] 


E(K,T)  = 


(1  -  e  C1T) 


0 


CE(l-e-ClT)  a2LE 
C1  C1C2 


[1  - 


C2e 


-C1T  n  _  2T 


+ 


c^e 


a2  al  a2  al 


—  [1  -  e  C2T] 


Q  +  Qp 

al  +  ^+UE(K) 

m 


Q  +  Q, 


a2  +  V 


m 


+  Ue(K) 
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2 . 3  Controllability  of  Discrete-Time  Models 

This  Chapter  will  be  concluded  with  a  discussion 
of  the  controllability  of  the  models  derived  earlier. 

THEOREM  [12] 

Given  a  continuous  system 

X  =  A  X  +  B  m 

where 

X  is  an  n  dimensional  state  vector 
m  is  an  r  dimensional  control  vector 
A  is  an  nxn  system  matrix 
B  is  an  nxr  control  matrix 

If  A  and  B  are  time  invariant  matrices  then  the 
system  is  completely  controllable  if  and  only  if  the 
controllability  matrix 

P  =  [B  |  AB  |  A2B  |  .  |  An_1] 

has  n  linearly  independent  column  vectors,  or  alternatively, 
if  the  matrix  P  has  a  rank  n. 

The  concept  of  controllability  for  discrete-time 
systems  is  similar  to  that  for  continuous  systems  except 
that  the  state  equations  are  difference  equations. 

The  condition  for  complete  controllability  of  a 
digital  system  is  such  that  given  the  system  equation 


X  (K  +  1 )  T 


A  (T)  X  (KT)  +  B  (T)  U  (KT) 


. 
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The  matrix 

[sq(T)  |  SX(T)  |  S2(T)  |  .  Sn_1(T)] 

Where 

s  (T)  =  A  (-T)  B  (T) 

S1(T)  =  A ( -2T) B (T) 


Sn-1(T)  =  A(“nT)B(T) 

and  T  is  the  sample  period. 

have  n  linearly  independent  column  vectors.  In  other  words, 

given  an  initial  state  X  there  exists  a  set  of  controls 

— o 

U(KT)  K  =  0,  1,  2  .  n-1  that  will  transfer  the 

system  to  any  desired  final  state  in  a  given  finite  number 
of  sample  periods. 

A  look  at  the  $ (T)  and  E (T)  matrices  of  the 
system  model  for  the  aeration  control  case  immediately 
reveals  that  the  system  is  not  completely  controllable, 
since  the  control  (KT)  has  no  effect  on  the  state  variable 
X2. 

This  is  not  a  set  back  since  the  main  objective 
is  to  control  the  concentration  of  dissolved  oxygen  in  the 
river.  Therefore  the  X2 (BOD)  variable  will  be  treated  as 
a  disturbance  term  in  the  design  of  a  controller  for  the 
system.  Only  the  controllable  part  of  the  model  will  be 
used  in  the  design  of  the  controller. 


■ 
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For  the  aeration  plus  effluent  discharge  case 

however,  the  $(K,T)  matrix  varies  between  sample  intervals 

since  it  is  a  function  of  U_ (K) ,  the  effluent  discharge 

control.  However,  the  system  is  controllable  since  a  set 

of  controls  [U_  (K)  U„ (K) ]  can  be  found  that  will  transfer 

a  r. 

the  system  from  any  initial  state  to  any  desired  final 
state  in  any  desired  number  of  sample  periods. 

It  will  also  be  noticed  that  the  controllability 
matrix  for  discrete  time  systems  is  a  function  of  the 
sample  period.  Thus  the  assertion  made  above  about  the 
controllability  of  the  system  is  on  the  assumption  that  a 
suitable  sample  period  can  be  found  such  that  the  linear 
independence  of  the  column  vectors  of  the  controllability 
matrix  is  not  violated.  The  choosing  of  the  sample  period 
for  the  systems  will  be  discussed  in  the  next  Chapter. 
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CHAPTER  3 


DESIGN  OF  DIGITAL  CONTROLLERS  FOR 
WATER  QUALITY 


3 . 1  Introduction 

In  this  chapter,  the  design  of  digital  controllers 
to  maintain  the  concentration  of  dissolved  oxygen  at  a 
prescribed  level  using  pole  assignment  technique  with 
output  feedback  is  discussed.  The  discrete  time  models 
derived  in  chapter  2  are  used  in  the  design. 

Consider  the  water  quality  control  system  shown 
in  figure  3.1  where  the  process  to  be  controlled  is  given  by 


X  (K+l )  T  =  $(T)x(KT)  +  E  (T)  U  (KT)  +  D  (T)  +  B(T)w(KT) 

(3-1) 


where 

x  is  the  state  variable  vector 
$  is  the  system  matrix 
U  is  the  control  vector 
E  is  the  control  matrix 

D (T)  and  B(T)w(KT)  are  disturbance  terms. 
K  =  0,1,2  .... 
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U  (KT) 

Controlled  Process 

x  (KT) 

N 

Figure  3.1  Water  Quality  Control  System 

Since  U(KT)  =  Nx(KT) ,  neglecting  the  disturbance 
terms,  equations  (3-1)  can  be  written  as 

x(K+l)T  =  $(T)x(KT)  +  E(T)Nx(KT) 

or 

x(K+l)T  =  (<f>(T)  +  E(T)N)  x  (KT)  (3-2) 

Taking  z  transforms  of  equation  (3-2) ,  we  have 
zX ( z )  -  x ( 0)  =  (<£>(T)  +  E (T) N) X ( z )  (3-3) 

Therefore 

X(z)  =  [zl  =  ($(T)  +  E(T)N)]"1  x  ( 0)  (3-4) 

where  I  is  the  unit  matrix. 

The  poles  of  the  closed  loop  system  transfer  matrix 
are  given  by  the  roots  of 

det [zl  =  ( $ (T)  +  E(T)N)]  =  0  (3-5) 


The  matrix  N  to  locate  these  roots  at  any  desired 
location  within  the  unit  circle  in  the  z  plane,  can  be 
determined  by  using  any  standard  method  of  pole  assignment 
[18,  19,  25]  . 
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3 . 2  Sampling  Interval 

As  pointed  out  in  chapter  2,  since  $  and  E  in 
equation  (3-5)  are  functions  of  the  sampling  period  of  the 
discrete  time  system,  it  is  necessary  to  determine  an 
appropriate  sample  interval  to  ensure  controllability  of  the 
systems.  For  this  purpose  the  sampling  theorem  will  be 
invoked  at  this  point. 

Sampling  theorem  [12] 

If  a  continuous  signal  contains  no  frequency  higher 

than  03  radians  per  unit  time,  the  signal  is  completely 
c 

characterized  by  the  values  of  the  signal  measured  at 
instants  of  time  separated  by 

T  =  7T  /  03  (3-6) 

c 

Case  1:  Artificial  aeration  control 

Using  the  values  used  by  Ramar  and  Gourishankar 
[18,  19]  viz 


al  = 

0.2 

a2 

0.32 

Q 

8.4  x  104 

qe 

2.8  x  104 

LE  “ 

20.2 

CE 

2.0 

► 

LA  “ 

0.0 

cs  = 

11.0 

db  = 

1.0 

Vra  = 

15.1  x  104  ) 

replacing  t  -  t  by  T 

in  equations 

(2-16),  (2-17), 

(2-18)  and  (2-19)  we  have 
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$  (T) 


-0. 9417T 
e 


0 


o  £-7  r  ~0 . 9417T  -1.0617T, 

2.67  Le  -e  ] 

-1. 0617T 
e 


(3-8) 


B  (T)  = 


0.59  (1-e  °*9417T) 


-0.178  (1  -  8.85e  °*9417T  + 


7 . 8  5e 


-1.0617T 


0 


0.52  (1-e  1'0617T) 


(3-9) 


D  (T)  = 


1.66 (1-e  °*9417T)  -  1.187  (1-8. 85e  °*9417T  + 
7.85e*1*0617T) 


3.49 (1-e  1-0617T) 


(3-10) 


E  (T)  = 


1.062  (1-e  °-9417T  ) 


0 


(3-11) 


Noting  that  the  open  loop  time  constant  of  the 

variable  x^  in  equation  (3-8)  is  ^ ^  or  1.0619  days  and 

that  the  open  loop  system  is  overdamped,  it  is  reasonable  to 

specify  the  highest  frequency  contained  in  the  continuous 

system  response  of  the  variable  to  be  4.0  radians  per  day. 

Therefore,  substituting  a)  =4.0  radians  per  day  into  equation 

o 

(3-6) ,  the  sampling  period  is 


3.1416 

4.0 


0.79  day 


This  sample  period  can  be  used  but  a  T=0.5  day 
is  more  convenient  (also  more  accurate).  Substituting  T=0.5 


. 
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into  equations  (3-8), 


<S> 


0.6245 

0 


(3-9) ,  (3-10)  and  (3-11)  we  have 

-0.0901 

0.5881 


(3-12) 


B 


0.2215 

0 


-0.0160 

0.2140 


(3-13) 


D 


0.5166 

1.4375 


(3-14) 


E 


.3987 

0 


(3-15) 


Substituting  for  $ ,  B,  D  and  E  from  equations  (3-12),  (3-13), 

(3-14)  and  (3-15),  equation  (2-21)  for  a  sampling  period  of 
T  =  0 . 5  becomes 


Xx (K+l ) 

.6245 

-0.0901 

(K) 

+ 

0.2215 

-0.0160 

Ci  (K) 

x2 (K+l) 

0 

0.5881 

x2  (K) 

0 

0.2140 

I-  (K) 

0.5166 

0.3987 

+ 

+ 

1.4375 

0 

(K) 


(3-16) 


. 
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Case  2 :  Artificial  aeration  plus  effluent  discharge 
as  controls 

It  was  shown  in  chapter  2  that  the  $  and  E  matrices 
for  this  case  are  dependent  on  the  effluent  discharge  control. 
Consequently  the  system  matrix  $  varies  from  one  sample 
period  to  the  next.  A  precise  computation  of  the  optimum 
sampling  period  would  be  difficult.  However,  since  we  can 
expect  the  response  of  the  variable  x-^  to  contain  higher 
harmonics  compared  to  the  single  control  case,  a  sampling 
period  less  than  0.5  day  would  seem  desirable.  For  convenience 
T=0.25  day  will  be  used.  The  suitability  of  the  choice  is 
confirmed  by  simulation  results  presented  later  in  this  chapter. 

Substituting  T=0.25  and  equation  (3-7)  into 
equation  (2-24)  and  noting  the  fact  that  the  nominal  value 
of  Q„;  Q  =  2.8  x  104  we  have 

hj  E 

x  (K+l )  =  $(K)x(K)  +  B  (K)  w  (K)  +  D  (K)  +  E(K)U(K)  (3-17) 


where 


<MK) 


0 . 7 902e  °-25ue(k)  -0 . 0075e  °*25UE(K) 


0 


0 . 7669e  °*25uE(K) 


(3-18) 


B  (K) 


(3-19) 


0.556  (1  -  0 . 7  902e 
0.9417  +  U 


-0 . 25Ue (K) 

^(K) 


where 
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-0.17  8  [1  -  (0.9731  +  0 . 9  4  2  U_  ( K )  )  e  °*25ue(k)] 

JCj 

(0.9417  +  U  (K) )  (1.0617  +  U„ (K) ) 


*21  =  ° 

0.556  [1  -  0.7669e"0,25UE (K) ] 
^22  1.0617  +  n  (K) 

£j 


D(K) 


1.571  [1  -  0.7902e  0*25ue(k)j 

0.9417  +  n (K) 

£j 

1.18  7 [1  -  (0.9731  +  0.194 2UE (K) )e-°*2^uE(K) 
(0.9417  +Ue(K))  (1.0617  +UE(K) 


3.71  (1  -  0 . 7669e  °*25ue(k) ) 


1.0617  +  U„  (K) 


(3-20) 


E(K) 


where 


a 


11 


a 


12 


1 — 1 
r—\ 

D 

al2 

°21 

a22 

[1  - 

0.7902e-°-25UE<K>] 

[0.9417  +  Ue(K) ] 

2(1  - 

■0.7902e-°-25UE'K) 

(0.9417  +  Ue(K) 


) 


(3-21) 


-0 . 25Ue (K) ] 


6 . 4  [1  -  (0.9731  +  Q.1942UE(K) )e 

[0.9417  +  Ue(K) ]  [1.0617  +  UE (K) ] 


a 


21 


0 
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a 


22 


2  0  [  1  -  0 . 7  669e  °*25uE(K) ] 

[1.0617  +  Ue(K) ] 


If  e  0*25UE(K)  expanded  in  a  Taylor  series  and 
the  first  two  terms  are  taken  (this  is  acceptable  if  UE(K) 
is  assumed  to  be  small  for  all  K  then 


:  °*25ue(K)  =  ].  -  0.25  U  (K) 

Ej 


(3-22) 


Substituting  equation  (3-22)  into  equation  (3-18) , 
(3-19),  (3-20)  and  (3-21)  and  assuming  that  th,  (K)  is  small 

£j 

compared  to  0.9417  for  all  K,  the  equations  become 


<MK) 


0.7902  -  0.197 6Ue (K) 
0 


-0.0075  +  0 . 0019Ue (K) 
0.7669  -  0 . 1917LL-,  (K) 

Ej 

(3-23) 


B  (K) 


D(K) 


0.1239  +  0.11671L  (K) 

Ej 

0 


.  0048  +  0 . 008  7Ut,  (K) 

III 

.1221  +  0 . 1004th,  (K) 

ill 

(3-24) 


0.3179  +  0.2711Ue(K) 
0.8148  +  0.6699Ue(K) 


(3-25) 


E  (K) 


0.2228  +  0.2098Ue(K) 
0 


0.2734  +  0 . lOlSU^ (K) 

ill 

4.3911  +  3 . 6112U  (K) 

E 


(3-26) 
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Equation  (3—17)  can  therefore  be  written  as 


x (K+l )  =  ($  +  $(K))  x (K)  +  [B  +  B (K) ]  w(K)  + 


[D  +  D(K)  ]  +  [E  +  E  (K)  ]  U  (K) 


where 


$  = 


0.7902 


0 


-0.0075 


0.7669 


$(K)  = 


-0.197  6Ue (K)  0.0019Ue (K) 

0  -0 . 1917Ue (K) 


B  = 


0.1239 


0 


0.0048 


0.1221 


B(K)  = 


0.1167Ue(K)  0.0087Ue(K) 

0  0.1004Ue(K) 


D  = 


E  = 


0.3179 

0.8148 

0.2228 

0 


D(K)  = 


0.2711Ue(K) 
0. 6699Ue (K) 


0.2734 


4.3911 


E  = 


0 . 2  098Ue (K)  0.1015Ue(K) 

0  3 . 6112Ue (K) 


(3- 


Equation  (3-27)  is  the  discrete  time  model  for  the  artificial 
aeration  plus  effluent  discharge  controls  case  for  a  sampling 
period  of  0.25  day. 


27) 


- 
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3.3  Digital  Controller  for  Artificial  Aeration  as 

the  Only  Control 


We  shall  now  design  the  controller  N  for  the  single 
control  case  which  will  take  care  of  deviations  from  the 
nominal  steady  state  values  of  the  state  variables.  These 
values  are 


X1 

=  6.0 

x2 

=  6.61 

Ci 

=  10.0 

L. 

=  6.0 

1 

°A 

=  0.55 

It  was  mentioned  earlier  that  the  input  aeration 
has  no  effect  on  the  BOD  output  of  the  system.  Now  the 
controllable  part  of  the  model  (equation  3-16)  is 

x1(K+l)  =  0.6245x1  (K)  -  0.0972x2  (K)  +  0.2215COOO 

-  0 . 016L . (K)  +  0.5116  +  0 . 3987LL  (K)  (3-29) 

1  r\ 


If  the  desired  output  DO  is  the  nominal  steady  state  value  of 
x1 ;  x  ,  the  error  between  the  desired  output  and  the  actual 
output  of  the  open  loop  system  is 


E1(K+1)  =  0.6245E1(K)  -  0.0972E2(K)  +  0.22156(^00 

-  0.01661^00  +  0.39876Ua(K)  (3-30) 


where 


* 
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E2 

x2 

X2 

«Ci 

=: 

C  . 

l 

-  C. 

l 

6L . 

1 

— 

L. 

l 

-  L . 
l 

5UA 

= 

UA 

,  < 

ID 

1 

In  order 

to 

obtain  a 

state  error, 

a 

third 

equation 

is 

E3 (K+l) 

= 

e3(] 

( 3-3  Oa ) 


(3-31) 


This  is  equivalent  to  the  integral  control  for  continuous 
systems . 

The  river  system  is  now  represented  by  the  augmented 

system. 


Ex (K+l) 

0.6245 

0 

E3 (K+l) 

-1 

1 

0.2215 

-( 

Ei 

(K) 

0.0972 

_E3 

(K) 

0 

E2(K)  + 


0 


-0.016 
0 


6C. (K) 

l 

6L . (K) 

l 


0.3987 

0 


6Ua(K)  (3-32) 


Let  6 U A  ( K )  =  N1E]_(K)  +  (K)  (3-33) 

Substituting  equation  (3-33)  into  equation  (3-32)  we  have 


Ex (K+l) 

(0.6245  +  0.3987N1) 

0.3987N2 

(K) 

E3 (K+l) 

-1 

1 

e3  (K, 

34. 


0.0972 

0.2215 

-0.016 

6C. (K) 

e2 (K)  + 

l 

0 

0 

0 

_6  L  i  ( K )_ 

The  characteristic  equation  of  this  system  is  given  by 


det 


(0.6245  +  0.3987^) 


1 


-0.3987N2 
(z  -  1) 


0 


that  is 

z2  -  (1.624  5  +  0.3 987N1) z  +  (0.624  5  +  0.3  987N1  +  0.3  98  7N2)  =  0 

(3-34) 


Noting  that  the  open  loop  system  is  overdamped, 
let  the  closed  loop  poles  of  the  system  represented  by 
equation  (3-33)  be  located  within  the  unit  circle  in  the  z 
plane  at 


z  =  0.3094  +  0 . 5488 j 

z  2  =  0.3094  =  0 . 5488 j 


(3-35) 


This  choice  of  closed  loop  poles  was  arrived  at  after  several 
computer  simulation  runs  with  different  pole  locations. 

The  characteristic  equation  of  the  closed  loop  system  for  this 
pole  location  is  therefore 


z2  -  0.6188z +  0.3969  =  0  (3-36) 

Equations  (3-34)  and  (3-36)  both  represent  the  characteristic 
equation  of  the  closed  loop  system,  therefore  comparing 
coefficients,  N-^  and  N2  are  found  to  be 


35. 


Nx  =  -2.5224 

(3-37) 

N2  =  1.9516 

Substituting  and  N2  into  equation  (3-33)  the  incremental 
control  law  is 

6Ua(K)  =  -2.5224E1(K)  +1.9516E3(K)  (3-38) 

3 . 4  Implementation  of  the  Artificial  Aeration  Control  Scheme 

and  Simulation  Results 

The  practical  implementation  of  this  control  law 
should  not  permit  negative  values  of  the  control  function, 

U  (K)  =  6U  (K)  +  U  ,  as  this  would  mean  removal  of  dissolved 

r\.  x\  xx 

oxygen  from  the  river.  It  therefore  follows  that  when  such 
conditions  occur,  that  is  U  (K)  is  negative,  the  control 

XX 

action  of  the  controller  is  not  desired  and  the  controller 
should  be  turned  off  until  the  DO  level  drops  down  to  the 
desired  level. 

This  constraint  means  that  the  control  function 
can  only  take  on  positive  values  (including  zero)  and  should 
be  set  equal  to  zero  if  the  function  is  negative. 

This  objective  can  be  achieved  by  inserting  a 
limiter  between  the  controller  and  the  zero  order  hold  in 
the  schematic  diagram  of  figure  2.1.  This  arrangement  is 
shown  in  figure  3.2.  This  scheme  was  simulated  on  a  digital 
computer  and  the  simulation  results  are  shown  in  figure  3.6 (a, 
b,  and  c)  for  an  assumed  input  DO  and  BOD  profiles  shown 
in  figure  3.3  (a  and  b)  .  The  DO  and  BOD  profiles  are  the 


J 
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staircase  approximation  of  those  used  by  Ramar  and 
Gourishankar  [18,  19]. 


Digital  Controller 


Figure  3.2  Artificial  Aeration  Control  Scheme 

3.5  Digital  Controller  for  Artificial  Aeration  Plus  Effluent 
Discharge  as  Controls 

The  design  of  a  controller  for  the  two  controls 
case  is  not  as  straight  forward  as  for  the  single  control 
case  since  the  state  transition  matrix  $  is  a  function  of  the 
effluent  discharge  control.  The  strategy  we  will  use  is  to 
split  the  model  represented  by  equation  (3-27)  into  a  fixed 
part  and  a  variable  part  and  design  the  controller  using  the 
fixed  part  of  the  model.  The  effects  of  the  variable  part 
of  the  model  can  then  be  minimized  by  carefully  choosing  the 
gain  matrix  N  of  the  controller  using  the  same  method  used 
in  the  last  section. 

Using  this  strategy,  equation  (3-27)  can  be 
written  as  the  sum  of  two  equations. 


' 


BOD  INPUT  IN  MG/L 
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o 

o 


(a) 


o 

o 


(b) 

Figure  3.3  Assumed  Upstream  DO  and  BOD  Profiles 


(a) 

DO 

(b) 

BOD 

OUTPUT  00  IN  tIG/L 

.00  2.00  4.00  6.00  0.00  10-00  12.00 
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o 

o 


(c) 


Figure  3.4 


Discrete-Time  System  Response  (Constant- 
Parameters)  Aeration  Control  Case 


(a)  Output  DO 

(b)  Output  BOD 

(c)  Artificial  Aeration 
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xf(K+l)  =  4>x(K)  +  EU  (K)  +  Bw  (K)  +  D  (3-39) 

~  ~  ~ 

x (K+l )  =  $x(K)  +  EU (K)  +  Bw (K)  +  D  (3-40) 


Equation  (3-39)  is  the  fixed  part  of  the  model  represented 
by  equation  (3-27)  and  equation  (3-40)  is  the  variable 
part . 


Using  the  pole  assignment  technique  and  considering 
only  equation  (3-39)  the  controller  is  designed  as  follows: 

If  we  once  again  assume  that  the  desired  output 
DO  is  the  nominal  value  of  x^;  x^ ,  then  the  error  between 
the  desired  output  x-^  and  the  actual  output  is 


E1 (K+l) 


E2 (K+l) 


0.7902 


0 


-0.0075 


0.7669 


E1(K) 


E2  (K) 


+ 


0.2228 


0 


0.2734 


4.3911 


SUa(K) 


VK) 


+ 


0.1239 


0 


0.0048 


0.1221 


6C  .  (K) 

l 


6L . (K) 

l 


(3-41) 


Where  E-.  ,  E0,  6C .  ,  SL.  and  SUa  are  defined  by  equation  (3-30a)  . 
1  2  l  l  A 


In  order  to  maintain  an  acceptable  steady  state 
error ,  a  third  variable  E3,  defined  by  equation  (3-31)  is 
introduced.  The  river  system  is  therefore  represented 
by  the  augmented  system 
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Ex  (K+l) 


e2  (K+l) 


E3 (K+l) 


0.7902 


0 


-1 


-0.0075 


0.7669 


0 


0 


0 


E1(K) 


e2  (K) 


e3(k) 


+ 


0.2228 


0 


0 


0.2732 


4.3911 


0 


sua(k) 


ue(k) 


+ 


0.1239 


0 


0 


0.0048 


0.1221 


0 


6C. (K) 

l 


6Li(K) 


(3-42) 


Let 


U(K)  = 


N. 


N. 


N. 


N 


E1  (K) 


e3  (K) 


(3-43) 


Substituting  for  U (K)  in  equation  (3-42)  we  have 


E1(K+1) 

(0.7902  +  0.2228N1 

-0.0075  0.2228N2+ 

Ej.  (K) 

E2 (K+l) 

= 

+  0.2734N3) 

4.3911N3 

0.2734N4 

0.7669  4.3911N4 

E2  (K) 

E3 (K+l) 

-1 

0  1 

3  <K)_ 

0.1239 


0.0048 


0 

0 


0.1221 

0 


6w  (K) 


(3-44) 


The  characteristic  equation  is  therefore  given  by 

z3  -  (2.5571  +  0.2  228N1  +  0.27  34N3) z2  + 

(2.1631  +  0.3937N,  +  0.2  2  28No  +  0.5160N-  +  0.2  7  34NJ  z  - 

1  2  3  4 

(0.606 +  0.1709N1 +  0.1709N2  + 0.2426N3 +  0.2426N4)  =  0 

(3-45) 


The  closed  loop  poles  of  the  discrete-time  system 
are  chosen  as 


z  =  0.7525 

(3-46) 

z 2  3  =  0.5207  +  0.5971 


In  making  the  above  choice  of  pole  location  we  have 
arbitrarily  chosen  the  real  pole  of  the  system  to  be  at  the 
point  S  =  -1.1  in  the  S  plane.  This  maps  into  the  point 
z  =  0.7525  in  the  z  plane.  The  complex  poles  were  then 
chosen  after  several  simulation  runs. 

The  characteristic  equation  for  this  pole  location 
can  therefore  be  written  as 


z  3  -  1.7939Z2  + 1.4114Z  -  0.4723  =  0  (3-47) 

Equation  (3-45)  and  (3-47)  are  the  same,  therefore  comparing 
coefficients  and  solving  the  resulting  equations  give 


I 


N][  =  -2.4401  +  5.2678N4 

N2  =  2.7977  -  0.5934N4  (3-48) 

N0  =  -0.8030  -  4 . 292  9N . 

3  4 

Since  N, ,  N0  and  N0  are  functions  of  N, ,  N.  should 
1  2  3  4  4 

be  chosen  such  that  the  variations  of  the  output  DO  due  to 
parameter  variations  are  as  small  as  possible. 

The  derivations  of  sensitivity  functions  for 
discrete-data  systems  are  available  in  the  literature  [10,  13, 
20] .  A  simple  method  for  investigating  the  sensitivity  of  the 
output  vector  of  a  multivariable  system  due  to  parameter 
variation  proposed  by  Soliman  and  Kevorkian  [20]  is  used  in 
finding  a  suitable  value  of  N4 . 

Soliman  and  Kevorkian  [20]  defined  the  sensitivity 
of  the  output  vector  of  a  multivariable  sample  data  system  to 
small  parameter  variations  as  the  ratio  of  the  norm  of  the 
output  vector  deviations  in  the  closed  loop  realization  of 
the  system  to  the  norm  of  the  output  deviations  in  the  open 
loop  realization  of  the  system.  Thus  the  sensitivity  of  each 
output  of  the  multivariable  system  to  small  parameter  variations 
can  be  written  as 

eCl (KT) 
e0l(KT> 

ec2(KT) | 

e02<KT>  I 


S  (KT) 
n. 


S  (KT) 
n~ 


. 
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S  (KT) 
n 

K 


ecK (KT) 

i^iKTy 


(3-49) 


Where  SRk(KT)  is  the  sensitivity  of  the  Kth  output  of  the 

multivariable  discrete  time  system. 

ec  (KT)  is  the  output  vector  deviation  due  to 
K 

parameter  variation  in  the  closed  loop  realization  of  the 

system  and  eQ  (KT)  is  the  output  vector  deviation  due  to 

K 

parameter  variation  in  the  open  loop  realization  of  the  system, 
Noting  that  a  good  feedback  design  will  ensure  that 


the  norm  of  the  output  vector  deviation  in  the  closed  loop 
realization  is  much  less  than  the  norm  of  the  output  vector 
deviation  in  the  open  loop  realization,  that  is 


eCl<KT) 


<  < 


e0l<KT) 


(3-50) 


ec2  <KT) 


<  < 


'O. 


(KT) 


etc . 

The  performance  of  the  two  controllers  in  a  closed  loop 
realization  of  a  system  for  the  same  controlled  process  could 
be  investigated  by  simply  comparing  the  norm  of  the  output 
deviation  due  to  parameter  variation.  The  objective  in  our 
case  is  to  find  the  value  of  that  will  minimize  the 
variations  of  the  output  DO. 
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A  graph  of  the  norm  of  the  deviation  of  the  output 
DO  from  the  desired  level  for  various  values  of  will 
indicate  the  value  of  for  which  the  norm  of  this  deviation 
is  minimum.  This  graph  is  shown  in  figure  3.5.  From  the 
graph  the  value  of  is  found  to  be  -0.1910. 

Substituting  for  in  equation  (3-48)  the  gain 
matrix  for  the  controller  becomes 

-3.4462  2.9110 

0.0169  -0.1910 


Therefore  the  control  law  is 

6U,  (K)  =  -3 . 44  62E-,  (K)  +  2.9110E~(K) 

A  1  J  (3-51) 

U  (K)  =  0.0169E1(K)  -  0.1910E3(K) 

From  u„(K),  the  quantity  of  effluent  discharged  into  the 
E 

river  is 


Q  (K)  =  15.1  x  104(Ue(K)  +  UE 


(3-52) 


For  reasons  stated  in  section  3.4,  U^OO  an<^ 

Q  (K)  can  not  take  on  negative  values,  therefore  the  control 
E 

law  would  have  to  be  nonlinear.  The  nonlinear  control 
law  is 


uA(K)  =  ^  >  0 
Ua(K)  =  o  ;  Y1  <  0 


Up.  (K)  =  Y2  ?  Y2  -  ° 


ue(r) 


y2  <  0 


(3-53) 


0 


-2^5  ~2D  IL5  "LO  I05 

VALUES  OF  N4  x  10 


Norm  of  the  Output  DO  Deviations  versus  N 


4 


Figure  3.5 


NORM  OF  OUTPUT  DO  DEVIATION  FROM  NOMINAL 
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where 

Y]_  =  -3.4462Ex(K)  +  2.9110E3(K)  +  U  (K) 

Y2  =  0.0169E1(K)  -  0.1910E3(K)  +  U  (K) 

The  quantity  of  effluent  discharge  is  therefore 

Qe  =  15.1  x  104  x  Y2  (3-54) 

Figure  3.6  is  a  schematic  diagram  of  the  closed 
loop  discrete  time  system  for  the  aeration  plus  effluent 
discharge  controls  case. 

Figure  3.7  (a,  b,  c,  and  d)  show  the  response  of 
the  system  for  input  DO  and  BOD  profiles  shown  in  figure 
3 . 3 (a  and  b) . 
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o 
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o 
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(c) 


(d) 


Figure  3.7  Discrete-Time  System  Response  (Constant- 

Parameters)  Aeration  Plus  Effluent  Discharge 
Controls  Case 


(a) 

Output  DO 

(b) 

Output  BOD 

(c) 

Artificial  Aeration 

(d) 

Effluent  Discharge 

CHAPTER  4 


SEASONAL  TEMPERATURE  VARIATION 

4 . 1  Introduction 

It  is  normally  necessary  to  specify  the  range  of 
weather  conditions  within  which  proper  operation  of  systems 
are  guaranteed.  For  example  a  look  at  the  specifications 
for  a  typical  electronic  calculator  shows  that  the  operating 
temperature  range  is  between  10°C  to  40°C  and  outside  this 
range  a  malfunction  is  to  be  expected.  It  is  therefore 
necessary  to  investigate  the  effect  of  seasonal  temperature 
variations  on  the  discrete-time  controllers  designed  in 
chapter  3.  An  investigation  of  the  effect  of  temperature 
on  the  continuous  time  systems  designed  by  Ramar  and 
Gourishankar  is  also  carried  out  in  this  chapter. 

4.2  Relationships  Between  System  Parameters  and  Temperature 

The  controllers  designed  for  maintaining  the 
level  of  DO  in  a  single  reach  of  a  river  system  have  been 
designed  on  the  assumption  that  the  system  parameters  a^ ,  a 2, 
o  ,  D  in  equation  (2-1)  are  all  constants.  In  1961,  a 
committee  on  sanitation  engineering  research  in  the  U.S.A.  [2] 
published  a  report  which  establishes  the  relationship  between 
water  temperature  and  stream  reaeration,  and  since  then  some 
researchers  have  been  taking  the  effect  of  temperature 
variation  on  the  river  system  parameters  into  account  in 
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their  analysis.  [3,  25]  Young  and  Beck  [25],  in  assessing 
the  suitability  of  their  model  suggested  that  a  term  dependent 
on  temperature  and  sunlight  conditions  be  added  so  that 
observed  data  could  be  correlated  with  values  obtained  by 
solving  the  system  equations.  Davidson  and  Bradshaw  [3]  have 
also  analyzed  the  effect  of  water  temperature  on  the  DO  and 
BOD  profiles  of  a  river  system  using  the  following  temperature 
dependent  parameters. 


(t) 

=  a±  (20°C)e°*  025t°c 

a2(T) 

=  a2  (20°C) e0, 0464t°C 

CS(T) 

=  Cs  (20°C)e"0,02lT°C 

VT) 

=  3  (25  -  0.028  ( t°C  -  30) 2) 

Where  3  is 

a  constant. 

t°C  is  the  temperature  of  the  water  in  degrees  Celsius 
a^(20°C)  is  the  reaeration  constant  at  20°C 
a2(20°C)  is  the  BOD  decay  rate  at  20°C 
Cg(20oC)  is  the  saturation  value  of  DO  at  20°C 
We  will  use  these  relationships  in  this  investigation. 

Since  data  for  the  water  temperature  for  which 
Young  and  Beck  [25]  calculated  the  values  for  the  various 
parameters  are  not  available,  a  seasonal  temperature  profile 
is  developed  using  the  1974  average  monthly  surface  temperature 
for  Great  Britain.  The  average  monthly  temperatures  for  the 
various  stations  and  the  overall  average  for  the  twelve 
months  of  1974  is  shown  in  Table  1. 


51 


U 

co 

I''- 

1 - 1 

CD 

H1 

rH 

00 

r- 

i — 1 

04 

in 

CM 

r-" 

cr\ 

CO 

c- 

CO 

CO 

W 

• 

Q 

H* 

in 

r- 

in 

CD 

in 

00 

CD 

00 

00 

r- 

r-" 

C" 

00 

i-" 

00 

00 

CT\ 

r- 

> 

IT) 

04 

CD 

m 

r- 

CO 

04 

in 

00 

in 

CO 

I — 1 

rH 

r- 

CO 

CD 

CM 

1 - 1 

CD 

O 

• 

2 

in 

in 

CD 

in 

in 

in 

m 

CD 

CD 

CD 

r- 

r" 

00 

04 

CD 

W 

S3 

Eh 

CD 

c- 

04 

CM 

in 

00 

m 

m 

CO 

in 

CD 

i — 1 

rH 

m 

CO 

1 - 1 

rH 

00 

Eh 

U 

• 

• 

o 

ID 

r- 

CD 

r" 

O' 

m 

r- 

04 

t"- 

r" 

r- 

00 

00 

r- 

00 

00 

04 

2 

H 

1 - 1 

i — 1 

CO 

CO 

i — 1 

o 

04 

r-" 

00 

04 

r- 

i — 1 

CM 

cre 

in 

CM 

r- 

CD 

CD 

Q 

04 

• 

W 

w 

O 

o 

l — 1 

o 

rH 

04 

I — 1 

o 

CM 

1 — 1 

i — 1 

CM 

CO 

CM 

CM 

CM 

CM 

CM 

i — 1 

Eh 

CO 

I— 1 

i — 1 

l — 1 

I — 1 

i — 1 

rH 

rH 

rH 

rH 

i — 1 

1 - 1 

rH 

rH 

l — 1 

1 - 1 

« — 1 

rH 

rH 

OS 

2 

O 

H 

04 

Eh 

o 

04 

in 

CO 

o 

04 

C0 

o 

rH 

in 

CM 

00 

CO 

00 

CM 

04 

H 

W 

o 

t 

Ph 

P 

D 

<N 

CM 

CO 

CO 

CM 

h1 

m 

in 

in 

CD 

in 

in 

in 

H1 

P 

C 

1 - 1 

1 - 1 

rH 

rH 

1 — 1 

I — 1 

rH 

i — 1 

• — 1 

rH 

I — 1 

I — 1 

rH 

l — 1 

rH 

rH 

rH 

1 - 1 

1 - 1 

CO 

D 

C 

PQ 

>H 

O 

O 

t"- 

o 

rH 

I — 1 

00 

in 

r- 

i — 1 

in 

l — 1 

CO 

00 

r" 

04 

CM 

Q 

01 

• 

H 

P 

D 

1 - 1 

CM 

CM 

CO 

■'T 

CM 

CO 

in 

in 

m 

CD 

in 

H1 

in 

H1 

< 

OS 

l"D 

1 - 1 

1 - 1 

1 - 1 

« — 1 

i — 1 

1 - 1 

rH 

1 — 1 

i — 1 

I — 1 

rH 

i — 1 

rH 

rH 

i — 1 

1 - 1 

rH 

1 - 1 

1 - 1 

&H 

O 

H 

£ 

OS 

W 

o 

1 - 1 

CD 

r" 

00 

r-~ 

o\ 

in 

CTi 

00 

in 

CO 

o 

m 

CM 

CO 

1 - 1 

04 

CQ 

w 

2 

• 

S3 

P 

o 

1 - 1 

rH 

rH 

CM 

o 

CM 

CM 

CO 

CO 

CO 

CO 

CO 

CO 

H1 

CM 

Eh 

Eh 

P 

1 — 1 

1 - 1 

i — l 

rH 

rH 

i — 1 

1 - 1 

rH 

1 - 1 

rH 

i — 1 

rH 

rH 

i — 1 

rH 

rH 

1 - 1 

t - 1 

1 - 1 

< 

OS 

1 — 1 

H 

Ol 

O 

1 — 1 

O' 

i — 1 

04 

["• 

04 

CO 

O 

CO 

r"- 

r- 

LO 

r-~ 

r" 

cn 

00 

H1 

w 

O 

Cm 

• 

P 

CD 

04 

o 

04 

O 

00 

O 

o 

1 - 1 

i H 

o 

o 

o 

rH 

o 

o 

o 

o 

O 

PQ 

os 

< 

!sj 

i — 1 

rH 

1 - 1 

rH 

rH 

rH 

i — 1 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

<C 

o 

Eh 

Eh 

P4 

c 

pcS 

Q 

in 

CM 

in 

in 

f" 

CO 

CM 

04 

in 

r — 1 

cr\ 

CTi 

CO 

r- 

in 

in 

c- 

w 

04 

• 

• 

os 

U 

< 

CD 

r- 

CD 

CD 

CD 

in 

CD 

00 

04 

00 

r- 

r-" 

CD 

00 

00 

00 

00 

04 

r- 

D 

H 

Eh 

Eh 

OS 

in 

r-'- 

00 

CO 

1 - 1 

CO 

in 

CD 

CM 

CO 

CO 

CO 

1 - 1 

r" 

00 

o 

m 

OS 

W 

P-i 

s 

w 

H 

P 1 

u 

£ 

in 

in 

in 

CO 

in 

in 

CD 

CD 

in 

in 

in 

CD 

m 

in 

CD 

r- 

in 

PQ 

CD 

CO 

rH 

r- 

in 

in 

CO 

00 

00 

in 

CO 

1 - 1 

[■" 

CO 

CO 

CM 

in 

Eh 

>H 

W 

• 

>h 

P 

S3 

P 

ffi 

EH 

2 

Cm 

in 

CD 

in 

CO 

in 

in 

CD 

in 

in 

in 

CO 

CD 

in 

CD 

CD 

r- 

uo 

2 

l — 1 

CM 

i-" 

o 

04 

CO 

00 

CD 

rH 

CD 

CO 

CJC 

CO 

o 

CO 

i — 1 

Eh 

O 

< 

• 

2 

S 

P 

in 

CD 

CD 

CD 

CO 

in 

in 

CD 

in 

in 

LO 

CD 

CD 

C" 

r- 

00 

CD 

O 

s 

H 

O 

<d 

U 

P 

S 

•H 

U 

a 

G 

P 

P 

g 

(D 

rd 

0 

G 

G 

G 

>4 

0 

CP  6 

-P  P 

P 

0 

rd 

0  P 

k 

2 

p 

rd 

a) 

U 

a) 

tr>  P 

03 

tn  tr>  P 

>.  p  a  p 

W 

o 

o 

a) 

G  rH 

G 

CO 

CD 

C 

G 

03 

G 

<D 

Sm 

(D 

G 

o 

H 

■  H 

0 

(DPP 

rd 

•H 

rd 

>1p 

•H 

•H 

(D 

0  rH 

0 

G 

0 

< 

Eh 

£ 

5m 

(D 

u 

G  P 

a  p 

<D 

U 

a  p  ^ 

P 

5 

a 

5m 

s 

a 

c 

5m 

0 

a) 

-H  M 

(D 

rH 

1 — 1 

g 

>H  P 

5m 

£ 

G 

rd 

rd 

G 

w 

Eh 

CL) 

4-> 

-H  P  P 

co 

CD 

a) 

<d 

fd 

•H 

rd 

O 

a) 

0 

H  rH 

0  rH 

> 

CO 

P 

CO 

Eh  <C  W 

W 

p 

PQ  >  2  CQ  is  O  2 

P 

o 

U 

PQ 

04 

<c 

52. 


J.C.  Ward  [23]  had  reported  that  seasonal  temperature 
variations  could  be  represented  by  a  sine  wave  of  the  form 

x  =  a  ( sin  (by  +  <J>)  )  +  7  (4-2) 


where 

t  is  the  temperature  in  °C 
a  is  a  constant 

b  =  360/365  or  0.986  degrees/day 

t  is  the  average  yearly  temperature  in  °C,  and 
y  is  the  number  of  days  starting  from  October  1.  That 
is,  y  =  l  for  October  1.  A  sine  wave  was  fitted  to  the  overall 
average  monthly  temperature  of  Table  1  using  the  least 
square  fit  method  (see  appendix  1).  The  values  of  a,  0  and 
t  are  obtained  as 

a  =  -4.3 

<J)  =  -23.1° 

t  =  9 . 1 9  °C 

The  equation  of  the  temperature  profile  is 

t  °C  =  -4 . 3 (sin ( 0 . 986y  -  23.1))  +  9.19°C  (4.3) 

Since  Young  and  Beck  collected  their  data  in  the 
summer  months,  it  is  reasonable  to  assume  that  the  values  of 
a  ,  a  ,  C  and  D  used  in  their  work  are  the  values  of  these 
temperature  dependent  parameters  evaluated  at  an  average 
summer  temperature.  Thus  if  we  assume  that  the  summer 
starts  on  the  first  day  of  June  and  ends  on  the  31st  day 

summer  temperature  can  be  calculated  thus 


of  August,  an  average 


■ 
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where 


0 


1 


e 


2 


sinede  +  9 . 19°C 


(by-j^  -  23.1)  radians 

(by2  -  23.1)  (y^-q)  radians 


and  are  the  number  of  days  between  October  1 
and  June  1  and  October  1  and  August  31,  respectively. 


=  244  days 

Y2  =  335  days 

0 ^  =  217.48°  or  3.7957  radians 

@2  =  307.21°  or  5.3618  radians 

_  -4.3  [-0.6047  -  0.7936] 

TSC  5.3618  -  3.7957 


t°C 


13 . 03  °C 


(4-4) 


Substituting  equation  (4-4)  into  equation  (4-1) 
and  using  a^(x)  =  0.2,  a2(T)  =  0*32,  Cs(t)  =  H*0  and 
D  (t  )=1 . 0  we  get 

D 


a1(20°C)  =  0.1444 

a2(20°C)  =  0.1748 

Cg ( 20  °C)  =  14.462 


(4-5) 


3 


0.0590 
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The  temperature  dependent  parameters  are  obtained  by  substituting 
equation  (4-5)  into  equation  (4-1) .  They  are 


a1  (t) 

=  0.1444e°*  025t°c 

a  2  ( x ) 

=  0.1748e°*0464T°C 

CS(T) 

=  14.462e_0*02lT°C 

VT) 

=  0.059  (25  -  0.028 

(t°C  -  30) 2 ) 


(4-6) 


These  temperature  dependent  parameters  are  only  valid  for  a 
temperature  range  of  between  5°C  and  35°C. 


4.3  Effect  of  Temperature  Variation  on  Discrete  Time  Systems 

Substituting  equation  (4-6)  for  the  system 
parameters  in  the  discrete  time  models  in  chapter  2,  the 
initial  conditions  for  each  quarter  of  a  year  are  calculated. 
These  initial  conditions  are  shown  in  Table  2. 

TABLE  2 

INITIAL  CONDITIONS  FOR  EACH  QUARTER 


Start  of  Each  Quarter 

x1  (0) 

x2<°) 

o 

ID 

o 

<  H 
Id 

October  First 

6.0 

6.68 

0.31 

0.0 

January  First 

6.11 

7.31 

0.0 

0.0 

April  First 

6.20 

7.12 

0.0 

0.0 

July  First 

6.0 

6.62 

0.70 

0.0 

The  discrete 

time  models 

are  then  simulated 

on  a  digital  computer 

for  a  period 

of  80 

days  (we 

have 

. 
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assumed  an  input  (upstream)  DO  and  BOD  profiles  for  an  80 
day  period  in  chapter  3)  in  each  quarter  of  a  year  using 
the  discrete  time  controller  designed  in  chapter  3  and  the 
initial  conditions  of  table  2.  The  results  obtained  are 
presented  in  figures  4.2,  4.3,  4.4,  4.5,  4.6,  4.7  and  4.8. 

A  discussion  of  this  results  will  be  found  in  chapter  5. 

4 . 4  Effect  of  Temperature  on  the  Continuous-Time  Systems 

After  substituting  equation  4-6  for  the  system 
parameters  in  the  continuous-time  model  of  the  river  system, 
the  continuous  system  was  also  simulated  on  a  digital 
computer**  for  a  period  of  80  days  in  each  quarter  of  the 
year.  The  continuous-time  controllers  designed  by  Ramar 
and  Gourishankar  [18,  19]  are  used  (see  appendix  2).  The 
initial  conditions  are  the  same  as  in  table  2.  The  results 
obtained  are  presented  in  figures  4.9,  4.10,  4.11,  4.12, 

4.13,  4.14  and  4.15.  A  discussion  of  these  results  and  a 
comparison  of  the  continuous-time  and  discrete-time  controller 
performance  is  given  in  chapter  5. 


**  A  continuous  system  modelling  program  (CSMP)  was  used  in  this 

simulation.  This  program  was  obtained  from  the  computing  science 
library . 
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(a)  (b) 


o 

o 


o 


(c)  (d) 


Figure  4.1  Seasonal  Temperature  Profiles  Used  for  the 

Investigation  of  the  Effect  of  Temperature  on 
the  Discrete  and  Continuous-Time  Systems 

(a)  First  Quarter  (Oct.  -  Dec.) 

(b)  Second  Quarter  (Jan.  -  March) 

(c)  Third  Quarter  (April  -  June) 

(d)  Fourth  Quarter  (July  -  Sept.) 


■ 

57 


(a)  (b) 


o 

o 


o 

o 


TIME.  DAYS 


(C) 


(d) 


Figure  4.2 


Discrete-Time  Response  (Temperature  Dependent 
Parameters) 


(a) 

Output 

DO 

(b) 

Output 

DO 

(c) 

Output 

DO 

(d) 

Output 

DO 

(Oct.  -  Dec.) 
(Jan.  -  March) 
(April  -  June) 
(July  -  Sept.) 
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o 

o 


®  o 


Figure  4.3  Discrete-Time  System  Response  (Temperature 

Dependent  Parameters)  Aeration  Control  Case 

(a)  Output  BOD  (Oct.  -  Dec.)  (b)  Output  BOD  (Jan.  -  March) 
(c)  Output  BOD  (April  -  June)  (d)  Output  BOD  (July  -  Sept.) 
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o 
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o 


Figure  4.4  Discrete-Time  System  Response  (Temperature 

Dependent  Parameters)  Aeration  Control  Case 

(a)  Artificial  Aeration  (Oct.  -  Dec.) 

(b)  Artificial  Aeration  (Jan.  -  March) 

(c)  Artificial  Aeration  (April  -  June) 

(d)  Artificial  Aeration  (July  -  Sept.) 
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Figure 


(a) 

(c) 


4  5  Discrete-Time  System  Response  (Temperature 


Dependent  Parameters) 
Discharge  Case 

Output  DO  (Oct.  -  Dec.) 
Output  DO  (Apr.  -  June) 


Aeration  Plus  Effluent 

(b)  Output  DO  (Jan.  - 
(d)  Output  DO  (July  - 


March) 
Sept . ) 
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(a) 


(c) 

Figure  4.6  Discrete-Time  System 

Dependent  Parameters) 
Discharge  Case 

(a)  Output  BOD  (Oct.  -  Dec.) 
(c)  Output  BOD  (Apr.  -  June) 


(b) 


o 


(d) 

Response  (Temperature 
Aeration  Plus  Effluent 


(b)  Output  BOD  (Jan.  -  Mar.) 
(d)  Output  BOD  (July  -  Sept.) 
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Figure  4.7  Discrete-Time  System  Response  (Temperature 

Dependent  Parameters)  Aeration  Plus  Effluent 
Discharge  Case 

(a)  Artificial  Aeration  (Oct.  -  Dec.) 

(b)  Artificial  Aeration  (Jan.  -  March) 

(c)  Artificial  Aeration  (April  -  June) 

(d)  Artificial  Aeration  (July  -  Sept.) 
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o 


(c)  (d) 


Figure  4.8  Discrete-Time  System  Response  (Temperature 

Dependent  Parameters)  Aeration  Plus  Effluent 
Discharge  Controls  Case 

(a)  Effluent  Discharge  (Oct.  -  Dec.) 

(b)  Effluent  Discharge  (Jan.  -  March) 

(c)  Effluent  Discharge  (April  -  June) 

(d)  Effluent  Discharge  (July  -  Sept.) 
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(a)  (b) 


o 

o 


o 

o 


TIME,  DAYS 


(c) 


(d) 


Figure  4.9  Continuous-Time  System  Response  (Temperature 

Dependent  Parameters)  Aeration  Control  Case 

(a)  Output  DO  (Oct.  -  Dec.)  (b)  Output  DO  (Jan.  -  Mar.) 

(c)  Output  DO  (Apr.  -  June)  (d)  Output  DO  (July  -  Sept.) 
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(c) 


(d) 


Figure  4.10  Continuous-Time  System  Response  (Temperature 

Dependent  Parameters)  Aeration  Control  Case 

(a)  Output  BOD  (Oct.  -  Dec.) 

(b)  Output  BOD  (Jan.  -  March) 

(c)  Output  BOD  (April  -  June) 

(d)  Output  BOD  (July  -  Sept.) 


66. 


a 

o 


o 

o 

c n 


o 

ID 

CM 


o 

__  o 
h-° 

cr~ 

01 

LU 

a: 

o 

UD 

o' 


O 

o 

^•00  20.00  40.00  "  50 .00  "  80.00 

TIME,  DAYS 


(a) 


(b) 


o 

o 


o 

o 


(d) 


Figure  4.11 


Continuous-Time  System  Response 
Dependent  Parameters) 


(a)  Artificial  Aeration 

(b)  Artificial  Aeration 

(c)  Artificial  Aeration 

(d)  Artificial  Aeration 


(Temperature 
Aeration  Control  Case 

(Oct.  -  Dec.) 

(Jan.  -  March) 

(April  -  June) 

(July  -  Sept.) 


OUTPUT  DO  IN  MG/L  OUTPUT  DO  IN  MG/L 

•  00  2.00  4.00  6.00  8.00  10.00  12.00  n  .00  2.00  4.00  6.00  8.00  10.00  12.00 
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Fiqure  4.12  Continuous-Time  System  Response  (Temperature 

Dependent  Parameters)  Aeration  Plus  Effluent 
Discharge  Case 


(a)  Output  DO  (Oct.  Dec.) 

(c)  Output  DO  (April  -  June) 


(b)  Output  DO  (Jan.  -  March) 
(d)  Output  DO  (July  -  Sept.) 
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(c)  (d) 

Figure  4.13  Continuous-Time  System  Response  (Temperature 

Dependent  Parameters)  Aeration  Plus  Effluent 
Discharge  Case 

(a)  Output  BOD  (Oct.  -  Dec.)  (b)  Output  BOD  (Jan.  -  March) 

(c)  Output  BOD  (April  -  June)  (d)  Output  BOD  (July  -  Sept.) 
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(a) 


(b) 


Figure  4.14  Continuous-Time  System  Response  (Temperature 

Dependent  Parameters)  Aeration  Plus  Effluent 
Discharge  Controls  Case 

(a)  Artificial  Aeration  (Oct.  -  Dec.) 

(b)  Artificial  Aeration  (Jan.  -  March) 

(c)  Artificial  Aeration  (April  -  June) 

(d)  Artificial  Aeration  (July  -  Sept.) 
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(a) 
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TIME,  DAYS 


(c) 


TIME.  DAYS 


(b) 
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(d) 


Figure  4.15  Continuous-Time  System  Response  (Temperature 

Dependent  Parameters)  Aeration  Plus  Effluent 
Discharge  Controls  Case 

(a)  Effluent  Discharge  (Oct.  -  Dec.) 

(b)  Effluent  Discharge  (Jan.  -  March) 

(c)  Effluent  Discharge  (April  -  June) 

(d)  Effluent  Discharge  (July  -  Sept.) 
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CHAPTER  5 


SUMMARY  OF  RESULTS  AND  CONCLUSIONS 


5 . 1  Summary  of  Results 

In  the  previous  two  chapters,  we  derived  discrete¬ 
time  controllers  for  maintaining  the  DO  level  of  a  polluted 
river  system  at  a  prescribed  level.  Two  cases  were  considered 
namely:  Artificial  aeration  control  only  and  Artificial 

aeration  plus  effluent  discharge  controls.  We  also  investi¬ 
gated  the  effect  of  seasonal  temperature  variation  on  these 
controllers  and  the  continuous-time  controllers  designed  by 
Ramar  and  Gourishankar  [18] .  The  results  are  summarized  below. 


1 .  Discrete-time  system  with  aeration  control  only 

(a)  System  model 

x  (K+l )  T  =  <E>  (T)  x  (KT)  +  B  (T)  w  (KT)  +  D  (T)  + 

E(T)Ua(KT) 

where  $ (T) ,  B (T) ,  D (T)  and  E (T)  are  defined  in  equations  (2-16), 
(2-17),  (2-18)  and  (2-19)  with  t  -  tQ  replaced  by  T. 


(b)  Control  function 


VK) 


Y  for  Y  >_  0 
0  for  Y  <  0 


where  Y  =  -2.5224  E± (K)  +  1.9516E3(K)  +  UA 
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2.  Discrete-time  system  with  aeration  and  effluent 
discharge  controls 

(a)  System  model 

X (K+l ) T  =  $(K,T)  x  (KT)  +  B(K,T)  w  (KT)  +  D(K,T)  + 

E(K,T)  U(KT) 

Where  $(K,T),  B(K,T),  D(K,T)  and  E(KfT)  are  defined 
in  equation  (2-14) . 


(b)  Controls 


VK) 


for  0 

0  for  Y^  <  0 


Where  Y1  =  -3.4462  E1 (K)  +  29110X3 (K)  +  UA 


qe(k) 


Y 2  for  Y 2  0 

0  for  Y 2  <  0 


where  Y9  =  15.1  x  104  [0 . 0169E,  (K)  -  0.1910X~(K)  +  U  ] 

£  1  j  E 


3 .  Effect  of  seasonal  temperature  variation 
(a)  Temperature  variation 


t  =  -4.3  [sin(0.986y  -  23.1)]  +  9.19°C 


(b)  Temperature  dependent  parameters 

The  parameters  a^,  a2,  C  ,  Dg  in  the  system 
models  are  replaced  by  temperature  dependent  parameters  for 
the  purposes  of  simulation.  These  temperature  dependent 
parameters  are  defined  in  equation  (4-6) .  In  the  design  of 
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the  controllers  however,  these  parameters  were  assumed  to  be 
constant  at  their  values  corresponding  to  the  mean  temperature 
during  the  summer  months  (June  -  August) . 

5 . 2  Discussion  of  Results 

(1)  One-control  scheme  ie.  aeration  control  only  - 

-  Effluent  discharge  kept  constant 

(a)  Constant  parameters 

The  results  obtained  for  the  aeration  control  case 
with  constant  parameters  are  shown  in  figure  3-4.  These 
results  are  similar  to  those  obtained  by  Ramar  and  Gourishankar 
[18]  for  the  continuous  controller. 

(b)  Effect  of  temperature  variation  on  the 
discrete  controller 

The  results  obtained  are  shown  in  figures  4.2, 

4.3  and  4.4.  The  results  show  that  the  amount  of  artificial 
aeration  required  to  maintain  the  desired  level  of  DO  at 

6  mg/£  decreases  as  the  temperature  of  the  river  decreases. 

This  is  due  to  the  increase  in  the  saturation  level  of  the 
dissolved  oxygen  in  the  river  and  the  increase  in  the  difference 
between  the  reaeration  rate  and  the  BOD  decay  rate.  These 
two  factors  produce  a  net  increase  in  the  amount  of  oxygen 
available  for  the  decomposition  of  the  waste.  Since  in  this 
case  the  effluent  discharge  is  kept  constant,  the  increase 
in  the  amount  of  oxygen  available  results  in  a  decrease  in 
the  amount  of  artificial  aeration  required. 
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(c)  Effect  of  temperature  variation  on  the 
continuous  controler 

The  results  obtained  for  the  continuous  controller 
are  shown  in  figures  4.9,  4.10,  and  4.11.  These  results 
show  unsatisfactory  DO  response,  especially  for  the  first 
three  quarters  (Oct.  -  June)  of  the  year.  This  is  because 
the  controller  was  not  turned  off  when  the  DO  level  is  higher 
than  the  desired  level.  A  non-linear  control  function  will 
be  required  in  order  to  improve  the  performance  of  this 
controller . 


(d)  Comparison  of  continuous-time  controller  and 
the  discrete  time  controller 

As  mentioned  in  section  (a)  above  the  results 
obtained  for  the  discrete-time  controller  with  constant  para¬ 
meters  compare  very  well  with  the  continuous-time  controller. 

The  number  of  observations  required  to  maintain  the  output 
DO  level  of  the  river  had  been  reduced  to  twice  per  day  without 
any  sacrifice  in  the  performance  of  the  system.  The  system 
could  easily  be  implemented  by  a  human  operator  using  a  set 
of  control  valves.  Since  the  effluent  discharge  is  kept 
constant  in  this  scheme,  it  might  be  necessary  to  store  the 
effluent  in  a  lagoon  so  as  to  smooth  out  the  variation  in  the 
production  of  effluent  from  the  sewage  plant  and  its  dumping 
into  the  river. 

2 .  Two-control  scheme  ie .  aeration  plus  effluent  discharge 

(a)  Constant  Parameters 

The  results  obtained  for  the  two  controls  case 


with  constant  parameters  are  shown  in  figure  3.7. 


These 


, 
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results  also  compare  very  well  with  those  obtained  by  Ramar 
and  Gourishankar  [19] .  While  the  volume  of  effluent 
discharged  for  the  continuous  system  varies  from  about  40,000 
to  0  cubic  meters  with  a  lot  of  fluctuations,  the  volume  of 
effluent  discharge  for  the  discrete  time  systems  never  goes 
below  20,000  cubic  meters  and  the  fluctuations  are  a  lot 
less.  However,  the  amount  of  artificial  aeration  required  to 
maintain  the  desired  level  of  DO  in  the  river  is  more  than  that 
required  in  the  continuous  case. 

(b)  Effect  of  temperature  variation  on  the 
discrete  time  controller 

The  results  are  shown  in  figures  4.5,  4.6,  4.7 
and  4.8.  These  results  show  that  the  amount  of  effluent 
discharged  into  the  river  increases  as  the  temperature  of  the 
river  decreases.  Thus  more  waste  will  be  assimilated  by  the 
river  during  the  winter  months.  The  reasons  are  as  stated 
in  section  5.2  subsection  1(b)  of  this  chapter.  This 
increase  in  the  amount  of  effluent  however,  is  also  accompanied 
by  an  increase  in  the  amount  of  artificial  aeration  required 
to  maintain  the  desired  level  of  DO  in  the  river. 

(c)  Effect  of  temperature  variations  on  the  continuous- 
time  controller 

The  results  for  this  part  of  the  investigation  are 
shown  in  figures  4.12,  4.13,  4.14  and  4.15.  The  results 
show  the  same  trend  as  for  the  discrete  time  controller  ie. 
the  amount  of  effluent  discharged  increases  as  temperature 
decreases,  however,  this  increase  in  effluent  discharge  is 
followed  by  a  decrease  in  the  amount  of  artificial  aeration 
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required  to  maintain  the  level  of  DO  at  6  mg/£ .  This  is 
different  from  what  was  obtained  for  the  discrete  time 
controller.  This  might  be  due  to  the  fact  that  the  sensitivity 
of  the  closed  loop  pole  to  parameter  variation  was  minimized 
in  the  design  of  the  continuous  controller,  this  was  not  done 
in  the  design  of  the  discrete  time  controller,  instead  the 
minimization  of  the  sensitivity  of  the  output  DO  to  parameter 
variation  was  used  in  determining  the  gain  of  the  discrete¬ 
time  controller. 

(d)  Comparison  of  the  discrete-time  system  and  the 
continuous  time  system 

The  results  show  the  same  trends  as  far  as  effluent 
discharge  is  concerned.  The  amount  of  artificial  aeration 
required  however,  show  opposite  trends,  but  it  should  also 
be  noted  that  the  fluctuations  in  the  amount  of  effluent 
discharged  is  greatly  reduced  for  the  discrete  time  system. 

The  amount  of  effluent  discharged  is  more  during  the  first 
half  of  the  80  day  period  than  the  second  half  for  the 
continuous  system. 

5 . 3  Conclusion 

Apart  from  the  increase  in  the  amount  of  artificial 
aeration  required,  it  had  been  shown  in  this  discussion  that 
the  performance  of  the  discrete-time  system  compares  very 
well  with  that  of  the  continuous  system.  With  the  reduced 
number  of  observations  required,  two  per  day  for  the  constant 
parameter  case,  and  four  per  day  for  the  temperature  dependent 
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parameter  case,  the  implementation  of  this  control  schemes 
are  much  simpler  than  the  implementation  of  the  continuous 
control  scheme.  As  stated  in  section  5.2  subsection  1(d) 
a  lagoon  might  be  required  for  temporary  storage  of  the 
effluent . 

5 . 4  Suggestion  for  Further  Research 

The  application  of  automatic  control  theory  to 
environmental  problems  is  still  in  the  development  stage. 

The  discussion  in  this  thesis  is  therefore  not  intended 
to  offer  a  solution  to  the  problem  of  organic  waste  pollution 
in  our  river  systems.  What  has  been  attempted  is  to  outline 
a  procedure  for  obtaining  a  solution  to  the  problem  using 
digital  control  techniques.  The  following  are  areas  for 
further  studies. 

(1)  In  chapter  two  digital  models  were  derived  from  a 
lumped  parameter  continuous  model,  since  the  validity  of 
the  parameters  of  these  models  is  restricted  to  a  limited 
temperature  range,  a  study  could  be  conducted  to  determine 
the  relationships  between  these  parameters  and  temperature 
for  lower  temperatures  up  to  and  including  the  freezing 
point  of  the  river,  followed  by  the  use  of  the  pole  assignment 
technique  to  regulate  the  output  variable  of  the  system  model. 

(2)  The  investigation  carried  out  in  chapter  4  had 
assumed  that  the  river  is  in  thermal  equilibrium  with  the 
atmosphere.  In  general, thermal  equilibrium  is  never  reached. 
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This  investigation  could  therefore  be  extended  to  include 
the  heat  exchange  mechanism  between  the  river  and  the 
atmosphere.  Also  the  effect  of  other  heat  sources,  such  as 
dumping  of  water  from  the  cooling  towers  of  power  generating 
plants  into  the  river  could  be  investigated. 
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APPENDIX  1 


DERIVATION  OF  LEAST  SQUARE  FIT  FOR  A  SINE  CURVE 

t  =  a  sin(bx  +  c)  +  t 

=  p  sinbx  +  q  cos  bx  +  t 

where 

c  =  tan-1 (q/p) 
a  =  /~p2  +  q2 

The  deviation  of  t  from  observed  data  t, 

d 

=  t  -  t, 
d 

=  p  sin  bx  +  q  cos  bx  +  t  -  t^  =  R 

The  sum  of  the  square  of  this  deviation  should  be  minimum 
for  least  square  fit. 

Thus 

ZR2  =  Z  [p  sin  bx  +  q  cos  bx  +  t  -  t^]  2  =  minimum 
differentiating  w.r.t.  P,  q  and  t  and  setting  the  derivatives 
to  zero  gives 

2Z  [p  sin  bx  +  q  cos  bx  t  t  -  t^]  sin  bx  =  0 

2Z  [p  sin  bx  +  q  cos  bx  +  t  -  t^jcos  bx  =  0 

2Z  [p  sin  bx  +  q  cos  bx  +  t  -  t^]  =  0 

This  is  equal  to 

PZ sin2  bx  +  q  Z  sin  bx  cos  bx  +  t  Z  sin  bx  -  E  t,  sin  bx  =  0 
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p  E  sin  bx  cox  bx  +  q  E  cos2  bx  +  t  E  cos  bx  -  E  cos  bx  =  0 

p  E  sin  bx  +  q  E  cos  bx  +  Et  -  Et^  =  0 

from  m  pairs  of  observed  data  for  x  and  the  last  equation 
can  be  written  as 

p  E  sin  bx  +  q  E  cos  bx  +  mt  -  Et^  =  0 

Therefore  solving  the  following  equations 
simultaneously  gives  p,  q  and  t. 

p  E  sin2  bx  +  q  E  sin  bx  cos  bx  +  t  E  sin  bx  -  E  t^  sin  bx  =  0 

p  E  sin  bx  cos  bx  +  q  E  cos2  bx  +  t  E  cos  bx  -  E  t^  cos  bx  =  0 

p  E  sin  bx  +  q  E  cos  bx  +  mt  -  E  t^  =  0 


APPENDIX  2 


CONTROL  FUNCTIONS  USED  IN  THE  SIMULATION  OF  THE 

CONTINUOUS  SYSTEMS 

Aeration  Control  Only 

UA(t)  =  -3.26  Ex(t)  +  8.82  E3  (t)  +  UA 

Aeration  Plus  Effluent  Discharge  Controls  Case 
U  (t)  =  -1.59  E1(t)  -  2.76  E3(t)  +  UA 
U  (t)  =  13.66  E]_(t)  -  52.11  E3(t)  +  UE 

The  controls  are  the  same  as  those  obtained  by  Gourishankar 
and  Ramar  [18,19]. 
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APPENDIX  3 


COMPUTER  PROGRAMS  USED  FOR  THE  SIMULATION  OF 


CLOSED  LOOP  SYSTEMS 

1 .  WPC 

(one-control,  constant  parameters) 

2.  WPC 2 

(2-control,  constant  parameters) 

3 .  SEW1 

(one-control,  temperature  dependent  parameters) 

4 .  SEW 

(2-controls,  temperature  dependent  parameters) 

5 .  WPT 

(continuous  system,  one-control,  temperature 

dependent  parameters) 

6 .  WPD 

(continuous  system,  two-controls,  temperature 

dependent  parameters) 

All  programs  are  written  in  Fortran.  WPT  and 

WPD  use  the  continuous  systems  modelling  program  (CSMP) 
available  in  the  computing  science  program  library.  All 
programs  were  run  on  the  Amdhal  470  V6  Computer  of  the 


University 

of  Alberta  Computing  Services  Department. 

A  listing  of  these  programs  is  included  in  this 

Appendix . 
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LIST  UPC 


>  1 

C  PROGRAM  TO  SIMULATE  CLOSED  LOOP  RIVER  SYSTEM ( CONSTANT  PARAMETERS) 

!:•  9 

REAL  XI (  162) 7  X2 ( 162 ) t DO  ( 162 ) -  BOD ( 162) >UA< 162) 7 TDO (80) 7 TBOD (SO) 

>  3 

REAL  TIM (162) 

4 

REAL  TIM ( 162 ) 

5 

T=0.5 

6 

TIM(1)=0.0 

>  7 

TIM ( 1 ) =0 . 0 

>  8 

TIM ( 161 ) =0.0 

>  9 

TIM  < 162 ) =20 . 0 

>  1 0 

XI  ( 16.1 )  =0 . 0 

>  1 1 

XI (162) =2.0 

>  1 2 

X2 ( 1 61 ) =0 . 0 

>  1 3 

X2( 162) =2.0 

>  1 4 

IJA  ( 161 )  =0 , 0 

>  1 5 

UA ( 162 ) =0 , 5 

>  1 6 

DO ( 161 ) =0.0 

>  1 7 

DO ( 162 )  =  1 . 0 

>  18 

BOD ( 161 ) =0.0 

>  19 

BOD( 162) =2.0 

>  20 

TIM ( 161 )  =0.0 

>  21 

TIM  < 162 ) =20 . 0 

9  9 

XI ( 161 )=0.0 

>  23 

XI ( 162) =2.0 

>  24 

X2 ( 161 ) =0 , 0 

*;:■  9  c; 

X2 ( 162 ) =2 . 0 

>  26 

UA( 161 )=0.0 

_■>  2  7 

UA( 162) =0.5 

>  28 

DO ( 161 ) =0 . 0 

>  29 

DO ( 1 62 ) =2 . 0 

>  30 

BOD ( 161 ) =0,0 

>  31 

BOD ( 1 62 ) =2 , 0 

>  32 

DO  30  K  = 1 .30 

>  33 

30 

READ ( 5  7  5000  )  TDO <  K ) y TBOD ( K ) 

>  34 

DO  40  K= 1 r  80 

>  35 

D0(2*K~1 >=TDO(K) 

>  36 

DO  (  2 * K  )  =  T D 0  (  K  )  +  (  TDO  (  K-f  1  )  -TDO  (  K )  )  /2 

>  37 

B0D(2*K-1)=TB0D(K) 

>  38 

40 

BOD ( 2*K ) =  T  B  0  D (  K  )  +  ( TBOD ( K  +  l ) -TBOD ( K ) ) /2 

>  39 

5000 

FORMAT (2F4. 1  ) 

>  40 

XX1=-0.9417*T 

>  41 

XX2=-1 , 0617*T 

>  42 

XI ( 1 )=6 

>  43 

UA ( 160 ) =9 , 0 

>  44 

XL=1 

>  45 

A=EXP ( XXI ) 

>  46 

TIM ( K+l ) =TIM  <  K ) +0 . 5 

>  47 

C=EXP ( XX2 ) 

>  48 

B=-2.67*(A-C) 

>  49 

D=0.6*< 1-A ) 

>  50 

E=-0.18*(1-8.85*A+7.85*C) 

>  51 

F=0.53*< 1-C) 

>  52 

0=1 .67* (1-A) -1 . 187* ( 1-8, 85*A+7. 85*0 

>  53 

H=3.5*( 1-C) 

>  54 

AA=(1-A)/0.9417 

5  5 

X2<  1  )  =  (F*6+H)./(  1-C) 

>  56 

U:L  =  (X1  ( 1  )*(  1-A)-B*X2(  1  )-D*10.0-E*6-G)/AA 

>  57 

UA ( 1 ) =U1 

>  58 

DO  50  K  = 1 7 1 5 9 

>  59 

T I M ( K  + 1 )=TIM(K)+0,5 

>  60 

UA (K+l ) =UA  <  K ) 

>  61 

ER-X1 (K)-Xl ( 1 ) 

>  62 

IF(ER.GT.O)  GO  TO  78 

63 

49 

X3=X3-ER 

>  64 

UA ( K ) =-2 . 522*ER+1 . 951 *X3+U1 

>  65 

IF ( UA ( K ) . LT . 0 )  UA<  K  >  =0 

>  66 

88 

GO  TO  66 

>  67 

78 

I F ( UA ( K ) » EQ , 0 )  GO  TO  88 

>  68 

GO  TO  49 

>  69 

66 

X2 ( K  +  l ) =C*X2 ( K ) +  F  *  B  0  D ( K ) +H 

>  70 

50 

XI (K+l )=A*X1 (K)+B*X2(K)+D*D0(K)+E*B0D(K)+G+AA*UA(K) 

>  71 

CALL  XPLOKTIM.Xl  »X2»UA»X1  .DOrBODfO.  1 ) 

>  72 

WRITE(6»6012) 

>  73 

6012 

FOR'MAT  (  '  INPUT  AERATION  CONTROL') 

>  74 

WRITE (6 7 200)  UA 

>  75 

WRITE ( 6  7  200 )  (XI (K) fK=l 7 160) 

>  76 

WRITE(6  7  6010 ) 

>  77 

6010 

FORMAT( 'OUTPUT  DO  MEASURED  TWICE  DAILY') 

>  78 

WRITE ( 6  7  200 )  (XI (K) 7K=1 7 160) 

>  79 

200 

FORMAT ( 20 ( 2X  7  F 4 . 2 ) ) 

>  80 

WRITE ( 6  7  601 1 ) 

>  81 

6011 

FORMAT ( 'OUTPUT  BOD  MEASURED  TWICE  DAILY') 

>  82 

WRITE (67 200)  (X2(K)  7K=:L  7 160) 

>  83 

WRITE ( 6  7  601 3 ) 

>  84 

6013 

FORMAT ( 'SYSTEM  PARAMETERS') 

>  85 

UR I TE ( 6  7  200 )  A  7  B  7  C  7  D  7  E  7  F  7  G  7  H 

>  86 

STOP 

>  87 

END 

O  / 

♦END  OF  FILE 


♦ 
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C  PROGRAM  TO  SIMULATE  CLOSED  LOOP  RIVER  SYSTEM  ( 2-CONTROL i  CONSTANT  PARAMETERS) 
REAL  X2( 162) ,  XI ( 162) > DO ( 160) >BOD( 160) >UA1 < 162) >UA2< 162) r TD0<80) »TB0D(30) 
REAL  TBOD ( 80 ) ?  TIM  ( 162 ) 

T  =  0 . 5 

TIM < 161 ) =0.0 
T IM  ( 162) =20 , 0 
XI < 161 >=0.0 
XI ( 162) =2.0 
X2< 161 ) =0.0 
X2( 162) =2.0 
UA.1  (161  )=().() 

UA1 (162) =0.5 
UA2 ( 161 ) =0 . 0 
UA2 ( 162 ) =1 . 0 
DO  30  K=1 y 80 

30  READ ( 5  ?  5000 )  TDO ( K ) ? TBOD ( K ) 

DO  40  K  =  1 >80 
D0(2*K-1)=TD0(K) 

DO ( 2*K ) =TDO (  K  )  +  ( TDO ( K+l ) -TDO ( K ) ) /2 
BOD ( 2KK-1 ) =TBQD ( K ) 

40  BOD ( 2*K ) =TBOD ( K ) + (TBOD ( K+l ) -TBOD ( K ) > /2 
5000  FORMAT  (  2F4  .  .1 ) 

TIM ( 1 ) =0 . 0 
Bl=— 3.44 
B2=2 .91 
B  3 = 0 . 0 1 6  9 
B4=-0, 1910 

WRITE ( 6  y  200 )  B1 ?  B2y  B3  p B4 

XI ( 1 )=6 

U2=0 

X3=0 

X4  =  0 

DO  50  K=1 1>  159 
TIM(K+1)=TIM  (10+0.5 
X  X 1  =  (  0 . 9  4 1 7  +  U  2  ) 

XX2=( 1 .0617+U2) 

A=EXP  (  -XXI  >X T  ) 

C=EXP ( -XX2&T ) 

B=-2 , 67.K  (  A-C  ) 

D=(0.56/XX1 )*< 1-A) 

E=-(  . 18/<XX1*XX2)  )*(  1-XX2*A/.  12+XX1*C/.  12) 

F=.56*(1-C)/XX2 

G  =  1 .57#(  1 -A >/XXl -.32*3.71/ (XXI *XX2>*(  1~XX2*A/.  12+XX1KC/.  12) 

H=3.71*( 1-C)/XX2 
AA= ( 1 -A ) /XX 1 

BB=2*( l-A) /XX 1-6.4/ <XX1*XX2)*(1-XX2*A/.12+XX1*C/. 12) 

CC=20*< 1-C)/XX2 

IF(K.GT.l)  GO  TO  45 

X2( 1 )=<F*6.0+H+CC*U2)/(1-C) 

U1=(X1 ( 1 )*(1-A)-B*X2< 1 )-D*10.0-E*6.0-G-BB*U2)/AA 
UA1 ( 1 ) =U1 
UA2( 1 )=. 19 

C  AERATION  CONTROL  INPUT 
45  UA1 (K+l )=UA1 (K) 

ER=X1 (K)-Xl (  1  ) 

IF(ER.GT.O)  GO  TO  73 
49  X3=X3-ER 

UA1 (K)=B1*ER+B2*X3+U1 
IF ( UA1 ( K ) . LT , 0 )  UA1 <  K ) =0 
GO  TO  250 

78  IF(UAKK)  .EQ.O)  GO  TO  250 
GO  TO  49 

C  EFFLUENT  DISCHARGE  CONTROL 

250  ER1=X1 (K)-Xl ( 1 ) 


tEND 


65 

UA2 (K+l )=UA2(K) 

66 

IF(ERl.LT.O)  GO  TO  202 

67 

201 

X4=X4-ER1 

68 

U2=B3*ER1+B4*X4 

69 

IF(U2.LT.-,19)  U2=-  .  1 9 

70 

UA2(K)=(0. 19+U2)*15. 1 

71 

GO  TO  350 

72 

202 

IF(UA2(K) .EQ.O)  GO  TO  350 

73 

GO  TO  201 

74 

350 

X2  (  K  +  l )  =C*X2  (  K  )  +  F  +  BOD  (  K  )  +H+CC5KU2 

75 

50 

XI ( K+l ) =A*X1 ( K ) +BKX2 ( K ) +D*BO ( K ) +E*BOD <  K ) +G  +  AA*UA1 <  K ) +BB*U2 

76 

UA1 (K+l )=UA1 (K) 

77 

UA2 ( K+ 1 ) =UA2 ( K ) 

78 

CALL  XPLOT (TIM»  XI >X2rUAl »UA2r 1 ) 

79 

WRITE ( 6  >  6010 ) 

80 

6010 

FORMAT( 'OUTPUT  DO  MEASURED  TWICE  DAILY') 

81 

WRITE ( 6 i 200 ) XI 

82 

200 

F0RMAT(10(2XvF8.2) ) 

83 

WRITE ( 6  >  601 1 ) 

84 

6011 

FORMAT ( 'OUTPUT  BOD  MEASURED  TWICE  DAILY') 

85 

WRITE (6 >200)  -  X2 

86 

WRITE  <  6  >  6012 ) 

87 

6012 

FORMAT  (  '  INPUT  AERATION  CONTROL.') 

88 

WR I TE ( 6  >  200 )  UA1 

89 

WRITE ( 6  r 6014 ) 

90 

6014 

FORMAT< 'EFFLUENT  DISCHARGE  CONTROL') 

91 

WR I TE ( 6  >  200 )  UA2 

92 

STOP 

93 

OF  FILE 

END 

* 
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C  PROGRAM  TO  SIMULATE  CLOSED  LOOP  RIVER  SYSTEM  ( 1 -CONTROL y  TEMPERATURE 
C  DEPENDENT  PARAMETERS) 

REAL  X2  ( 322 ) ,  XI  (322) »XDO< 160) rXBOD< 160) ylJAl (322) ylJA2(322> >TD0<80) 
REAL  TBOD ( 30 ) *  TEMP (322) -TIM (322) *  BOD (320) > DO (320) 

T=0 . 25 

TIM ( 322 ) =1 5 . 01 
XI (321 )=0.0 
XI (322>=2,0 
X  2  (  3  2 1 )  =  0 . 0 
X2 ( 322 ) =2 . 0 
UA1 ( 321 ) =0 . 0 
UA1 (322) =0.5 
UA2 ( 32 1 ) =0 . 0 
UA2 ( 322 ) -1 . 0 
TEMP ( 32 1 ) =0.0 
TEMP (322) =2. 5 
TIM ( 1 ) =0 . 0 
DO  30  K=lr80 

30  READ ( 5  f 5000  )  TOO ( K ) v TBOD (K) 

DO  40  K-1,80 
XD0(2*K-1)=TD0(K) 

XD0( 2*K ) =  TDG ( K ) + (TDO (K+l ) -TDO (K ) )/2 
XB0D(2*K-1)=TB0D(K) 

40  XBOD ( 2#K ) =  TBOD (  K )  +  ( TBOD ( K+ 1 ) -TBOD ( K ) ) /2 
DO  44  K=1 f 160 
D0(2*K-1 )=XDO(K) 

DO  (  2*K  )  =XDO  (  K  )  +.<  XDO  (  K+ 1 )  -XDO  (  K  )  )  /2 
BOD ( 2*K- 1 ) =XBQD ( K ) 

44  BOD ( 2*K ) =XBOD ( K ) + ( XBOD ( K+l ) -XBOD ( K ) ) /2 
50PO  FORMAT ( 2F4 . 1 ) 

Bl=-2.5221 
B2=l .9514 
B  3  =  0 . 0 
B  4  =  0 . 0 
DR =6 . 0 
U2=0 
X3=0 
X4=0 
X5=0 

READ (7 » 560)  XJyXl(l)rUl 
560  FORMAT ( 3F5 , 1 ) 

DO  50  K - 1 r  3 1 9 

T I M ( K+ 1 ) =T I M ( K ) +0 , 25 

XT=-4 . 30*SIN (0.01 72* ( X J ) -0 . 40 ) +9 . 1 9 

TEMP ( K ) =XT 

A 1=0. 1444* EXP (0.025* XT) 

A2=0 , 1748*EXP ( 0 . 0464*XT ) 

CS=1 4 . 462*EXP ( -0 . 0210*XT ) 

DB=0. 0590* (25. 0-0.023* (XT-30)* (XT-30) ) 

XX1=A1+0.7417+U2 
XX2=A2+0 . 7417+U2 
CD=XX2-XX1 
A  =  EXF'(-XX1*T) 

C=EXF'(-XX2*T) 

B=- ( A2/ ( A2-A1 ) ) * ( A - C ) 

D= ( 0 . 56/XX1 ) * ( 1 -A ) 


♦  END 


57  E=-0 , 56*A2* ( 1- ( XX2/ ( XX2-XX1 ) ) *A+ ( XXI/ ( XX2-XX1 ) )*C) 

58  F=(0.56/XX2)*( 1-C) 

59  G=(  (  Al*CS-DB  +  0 . 37  )  /XXI  )  *  ( 1  -  A  )  ~A2*3 . 71*  ( 1.  -  (  XX2/CD  )  *A+  (  XXI  /CD  )  *C  ) 


60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

31 

32 

33 

84 

85 

86 
87 


H=3.71*( 1-C)/XX2 
AA= ( 1 -A ) /XXI 

BB=2*( 1-A)/XX1~6.4/(XX1*XX2)* ( l-XX2*A/» 12+XX1 *C/ , 12) 

CC=20*( 1-C)/XX2 

IF(K.GT.l)  GO  TO  45 

X2( 1 )=<F*6.0+H+CC*U2)/(1-C) 

IF(Ul.LT.O.O)  U1=0,0 
UA.1  ( 1 )  =U  1 

C  AERATION  CONTROL  INPUT 
45  ER=X1 ( K ) -DR 

IF(ER.GT.O)  GO  TO  73 

49  X3=X3-ER 

UA.l  ( K)  =  B1*ER  +  B2*X3  +  U1 
IF(UA1(K) .LT.O)  UA1(K)=0 
GO  TO  350 

78  IF(K.GT.l)  GO  TO  88 

IF(UA1(K) .EQ.O)  GO  TO  350 
GO  TO  49 

88  IF (UA1 ( K  —  1 ) .EQ.O)  GO  TO  351 
GO  TO  49 

351  UA1 (K)=UA1 (K-l ) 

350  X2 ( K+ 1 ) =C*X2 ( K ) +F*BOD ( K ) +H+CC*U2 

XI ( K  +  l ) =A*X1 ( K  )  +B*X2 ( K ) +D*DO ( K ) +E*BOD ( K ) +G+AA*UA1 ( K ) +BB*U2 
WRITE <  6  » 200)  X 1 ( K ) y  X2 ( K ) «  UA 1 ( K ) *  X3  y  XT 
200  FORMAT ( 10(2XyF8,2)  ) 

50  X J=X J+0 , 25 
WRITE (6- 200)  XJ 
STOP 


88  END 

OF  FILE 


* 


88 
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C  PROGRAM  TO  SIMULATE  CLOSED  LOOP  RIVER  SYSTEM  ( 2- CONTROL ..  TEMPERATURE 
C  DEPENDENT  PARAMETERS) 

REAL  X2 ( 322 ) ,  X1C322) »XD0(160) yXBOD(160> »UA1(322) *UA2<322> » TDO (30) 
REAL  TBOD ( 80 ) , TEMP ( 322 >  y IIM ( 322 ) y BOD ( 320 ) » DO  ( 320 ) 

T~0 . 25 

T I M  (  3  2  2  )  =  1 5 . 0 1 
XI (321 )=0.0 
XI (  3  2  2  )  =  2 . 0 
X2 ( 321 ) =0 . 0 
X2 ( 322 )  =2 . 0 
UA1 (321 >=0.0 
UA1 ( 322  >=0.5 
IJA2(321  >=0,0 
UA2 ( 322  >  =1 , 0 
TEMP ( 321 >  =0 . 0 
TEMP (322) =2. 5 
TIM ( 1 >=0.0 
DO  30  K  = 1 ?  8 0 

30  READ ( 5 1 5000  )  TDO ( K > » TBOD ( K ) 

DO  40  K=l»80 
XDO ( 2+K-l >  =TDO ( K ) 

XDO  (  2  +  K  )  =TDO  (  K  >  +  (  TDO  (  K  +  1  >  -TDO  (  K  )  )  /2 
XBOD (2+K-l >  =TBOD ( K > 

40  XBOD < 2+K >  =TBOD ( K >  +  ( TBOD ( K  +  l > -TBOD  (  K ) >  /2 
DO  44  K= 1 y 1 60 
D0(2*K-1 >=XDO(K> 

DO  (  2+K  )  =XDO  <  K  >  +  <  XDO  (  K+ 1  >  -XDO  (  K  >  )  / 2 
BOD ( 2  +  K-l ) =XBOD ( K > 

44  BOD <2+K)=XB0D<K)+ (XBOD (K+l >-XBOD(K) )/2 
5000  FORMAT ( 2F4 . 1 ) 

B 1 =-3 . 44 
B2=2 .91 
B3=0 .0169 
B4=-0. 1910 
D  R  =  6 . 0 
U2=0 
X3=0 
X4=0 
X5  =  0 

READ ( 7  y  560  )  XJ y XI ( 1 > y U1 
560  FORMAT ( 3F5. 1 > 

DO  50  K= 1 y  31 9 

TIM (K+l )=TIM (K>+0,25 

XT=-4 . 30+S IN ( 0 . 0 1 72+  <  X  J  )  -0 . 40  >  +9 . 1 9 
TEMP ( K  >  =XT 

A1=0, 1444+EXP (0.025+ XT) 

A2=0. 1748*EXP(0.0464*XT) 

CS=1 4 . 462+EXP ( -0 . 02 1 O+XT ) 

DB=0 . 0590+ ( 25 . 0-0 . 028+ ( XT-30 )  +  ( XT-30 >  > 

XX1=A1+0.7417+U2 
XX2=A2+0 , 7417+U2 
CD=XX2-XX1 
A=EXP ( -XX1+T ) 

C=EXP ( -XX2+T > 

B=-( A2/( A2-A1 > )+( A-C) 

D=(0.56/XX1 >+( l-A) 

E=-0.56+A2+< 1-(XX2/(XX2-XX1 > ) +A+ ( XXI/ ( XX2-XX 1 ) >*C> 

F= ( 0 , 56/XX2 )  +  ( 1~C  > 

G=( ( Al+CS-DB  +  O  .37) /XXI >  +  ( l-A ) -A2  +  3 . 71  + ( 1 - ( XX2/CD > +A+ ( XX1/CD ) +C ) 
H=3,71  +  (  1-0/XX2 
AA= ( 1 -A ) /XX 1 

BB=2+( l-A>/XXl-6,4/(XXl+XX2)+( 1-XX2+A/ , 12+XX1+C/, 12) 

CC=20+( 1-C>/XX2 


' 


>  64 

IF(K.GT.l)  GO  TO  45 

>  65 

X2(  1  )=<F*6»0+H+CC*U2)/C1-C> 

66 

IF ( U1 , LT . 0 . 0 )  U1=0.0 

>  67 

UA1 ( 1 )=U1 

>  68 

UA2( 1 )=. 19 

>  69 

C 

AERATION  CONTROL  INPUT 

>  70 

45 

ER=X1 (K)-DR 

>  71 

U  A 1 ( K  T 1 ) -UA1 ( K ) 

>  7  ** 

IF(ER.GT.O)  GO  TO  78 

>  73 

49 

X3=X3-ER 

>  74 

UA1 <K)=B1*ER+B2#X3+IM 

>  75 

IF(UA1(K) *LT.O)  UA1(K)=0 

>  76 

GO  TO  250 

>  77 

78 

IF(K.GT.l)  GO  TO  88 

>  78 

IF (UA1 (K) *EQ.O)  GO  TO  250 

>  79 

GO  TO  49 

>  80 

88 

IFCUAi CK-1 > ,EQ.O)  GO  TO  250 

>  81 

GO  TO  49 

>  82 

C 

EFFLUENT  DISCHARGE  CONTROL 

>  83 

250 

ER1=X1  (K')-DR 

>  84 

UA2 ( K+l ) -UA2 (K> 

>  85 

IF ( ER1 . LT. 0  )  GO  TO  202 

>  86 

201 

X4=X4-ER 1 

>  87 

U2=B3*ER1+B4*X4 

>  88 

IF  <  U2 . LT  *-*19)  U2--.19 

>  89 

UA2 ( K )  =  ( 0 . 1 9+U2 ) *15  « 1 

>  90 

GO  TO  350 

>  9 1 

202 

IF(K.GT.l)  GO  TO  99 

>  92 

IF ( UA2  <  K ) , EQ . 0 )  GO  TO  350 

>  93 

GO  TO  201 

>  94 

99 

IF(UA2(K-1> .EQ.O)  GO  TO  350 

>  95 

GO  TO  201 

>  96 

350 

X2<K-H)=C*X2<K)+F#B0D<K)+H+CC#U2 

>  97 

XI  <K+1)=A*X1  <K)+B3KX2<K) +11*110  <K)+E#BOD(K)+G  + AA#UA1  <K)+BB*U 

>  98 

WRITE ( 6 i 200 )  X1(K) yX2(K) rUAl (K) r UA2(K) «XT 

>  99 

200 

FORMAT (10<2XrF8«2) ) 

>  100 

50 

X J=X J+0  *  25 

>  101 

WRITE  <  6  y  200 )  XJ 

>  102 

STOP 

>  103 

END 

#END  OF  FT  LI 

E 

w 


90 


LIST  WF'T 

1  LABEL  WATER  POLLUTION  CONTROL-  TEMPERATURE  EFFECT 

1.2  LABEL  SIMULATION  OF  CONTINUOUS  SYSTEM  (1-CONTROL) 

>  2  PARAMETER  X=1.0»  IC1=6.0?  UAB=0.31 

>  3  INITIAL 

>  3.1  TEMP=-4.30*SIN(0.0172*(X+TIME>-0.  -40)19. 19 

3.2  A1=0. 1444*EXP(0,025*TEMP) 

>  3.3  A2=0. 1748*EXF'(0.0464*TEMP) 

>  3.4  CS=14.462*EXP(-0,0210*TEMP) 

>  3 . 5  DB=0 . 0590* (25,0-0. 028* ( TEMP-30 ) *  <  TEMP-30 ) ) 

>  5  I C2= ( 0 . 556*6 , 0+3 , 7 ) / ( A2+0 . 74 ) 

>  6  I C  3  =  0 . 0 

>  S  DYNAMICS 

>  10  TEMF'=-4 , 30*SIN  (  0 . 0 1 72*  (  X  +  TIME  )  -0 . 40  )  +  9 . 1 9 

10.8  A1=0. 1444*EXF'  (0,025*  TEMP) 

>  11.6  A2=0,1748*EXP(0.0464*TEMP) 

>  12.4  CS=1 4 . 4<S2*EXP  (  -0 . 0210*TEMP  ) 

>  13,2  DB=0 . 0590* ( 25 . 0-0 . 028* ( TEMP-30 ) * ( TEMP-30 ) ) 

>  14  B0D=3  ,  *STEP  (  0  ,  )  +  ,4*RAMP<4.  )  -  .  4*RAMP  ( 14  .  )  + .  25*RAMP  (  24  .  )  .  .  , 

>  15  -1 ,25*RAMP(2S. ) +1 , 075*RAMP ( 31 . ) - . 075*RAMP ( 51 . )+. . . 

>  16  2 , 5*RAMF'  (  53  .  ) -5  ,  *RAMP  (  55  .  ) +2. 5*RAMP  <  57  .  )  -  .  2*RAMP  (  60  .  )  .  .  . 

>  17  + . 5*RAMP ( 70 .  ) 

IS  D0  =  10.*STEP(0.  )  FRAMPdO.  ) -RAMP  <  10 , 5  >-2  .  *RAMP  (21 .  )  +  .  .  . 

>  19  2  .  *RAMF'  (21.5)  +  RAMF'  (31,  )-RAMP(31 ,5  ) -3  .  *RAMP  (  41  .)+... 

20  3  .  *RAMF'  (41.5)  +RAMF'  (  51 .  )  -RAMP  (51.5) 

>  20.1  ERR-Xl-6 . 0 

>  20,2  X3-- INTGRL  ( IC3  r  -ER'R  ) 

>  20.3  X3B=-3 . 26*ERR 

>  20.4  UA1=3.82*X3+X3B+UAB 

>  20.42  UA=L..IMI  T  (0.0*  1 000  ♦  0  »  UA 1  ) 

>  20 . 5  X1D=- ( A 1+0 . 74 ) *Xl-A2*X2+0 . 556*D0+A1*CS-DB+ . 37+UA 

>  21  X2D=~ ( A2  +  0  »  74 ) *X2+0 . 556*B0D  +  3 , 7 

>  22  X1=INTGRL( IC1 rXlD) 

>  23  X2=INTGRL ( IC2  y X2D ) 

>  24  TIMER  DELT=0 . 1 r  F I NT  I M  =  80 t PRDEL=0 . 5 

>  25  PRINT  XI y  X2 >  TEMP  r  UA 

>  26  END 

>  26.1  PARAMETER  X=93,0»  181=6.61?  UAB=0.0 

>  26.2  END 

>  26.3  PARAMETER  X=183.0?  IC1=6.20?  UAB=0.0 

>  26,4  END 

>  26.5  PARAMETER  X=274,0?  IC1=6.0?  UAB=0.70 

>  26.6  END 

>  27  STOP 

>  28  ENDJOB 

0--END  OF  FILE 

* 


/> 
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LIST  WPD 

1  LABEL  WATER  POLLUTION  CONTROL-  TEMPERATURE  EFFECT 

>  1.2  LABEL  SIMULATION  OF  CONTINUOUS  SYSTEM  (2-CONTROL) 

>  2  PARAMETER  X=1.0r  ICl=6.0r  UAB--0.31 

>  3  INITIAL 

>  4  TEMF'=-4 . 30#SIN ( 0 . 0172# ( X+TIME ) -0  »  40  >  +  9 . 1? 

>  5  A 1  =  0 . 1 444*EXF'  (  0 . 025.YTEMF' ) 

>  6  A2  =  0  , 1748YEXF' (0.0464*  TEMP) 

>  7  CS=1 4 . 462*EXP ( -0 . 0210*TEMP ) 

>  8  DB=0 . 0590* ( 25 . 0-0 . 028* ( TEMP-30 ) * ( TEMP-30 ) ) 

>  10  IC2= ( 0 . 556*6 . 0+3 . 7 ) / ( A2+0 . 74 ) 

>  11  I  C3=0 . 0 

>  12.2  UEB=2 . 8 

>  12.4  0M  =  1 5 . 1 

>  13  DYNAMICS 

>  14  TEMP=-4.30*5IN(0,0172*(X+TIME>-0,40)+9. 19 

>  15  A  1=0. 1444*EXP(0.025*TEMF'> 

>  16  A2=0.1743*EXP(0.0464*TEMP> 

17  CS=14,462*EXP<~0,0210*TEMF') 

>  IS  DB=0.0590*(25,0-0.028*( TEMP-30) *( TEMP-30) ) 

>  1 9  B0D=3  .  *STEF'  (  0  .  )  +  .  4*RAMP  (  4  .  )  -  .  4# RAMP  (14.)  +  .  25:!<RAMP  (  24  ,  )  .  .  . 

>  20  -1.25*RAMP<28. )+l ,075*RAMP(31 . > - . 075*RAMP ( 51 .)+... 

21  2 . 5*RAMP  ( 53  .  > -5.  *RAMP(  55 .  ) +2 . 5*RAMP  (57 .  )  - .  2*R'AMF'  ( 60  .  )  .  .  . 

>  22  +.5*RAMP<70. ) 

>  23  D0=10  .  *STEP  (  0  ,  )+RAMP(10.  ) -RAMP  ( 10 . 5 ) -2  .  *RAMP  ( 21 .  )  +  .  .  . 

>  24  2.*RAMP(21 .5)+RAMP(31 . ) -RAMP ( 31 . 5 ) -3 . *RAMP ( 41 ,)+... 

>  25  3 ,  *RAMF'  (41,5)  +RAMP  (51 . )  -RAMP  (51.5) 

>  26  ERR=X 1-6.0 

>  27  X3=INTGRL(IC3r-ERR) 

>  29  UA1=-2.73*X3-1 ,59*ERR+UAB 

>  30  UA=LIMIT(O.Oy 1000.0 yUAl ) 

>  30.2  UE1=13.66*ERR~52. 11*X3+UEB 

>  30.4  UE  =  LIMIT (  0  ♦  0  r 1 ,0E12»UE1  ) 

>  31  XlD=-(Al+0.556+UE/VM)*Xl-A2*X2+0,556*D0+Al*CS-DB+2*UE/VM+UA 

>  32  X  2  D  =  - ( A2  +  0 . 556  +  UE/VM ) *X2+0 . 556*B0D+20*UE/UM 

>  33  X1=INTGRL(IC1»X1D) 

34  X2  =  INTGRL ( IC2  ?  X2D ) 

>  35  TIMER  BELT =0 . 1 »  F  INTIM=30  t  OUTDEL  =  .L  .  0  y  PRDEL=,5 

>  36  PRINT  X1*X2»UA»UE»X3*TEMP 

>  37  END 

>  37,1  PARAMETER  X=93.0y  IC 1=6. 61k  UAB=0.0 

>  37.2  END 

>  37.3  PARAMETER  X=133.0,  ICl=6.20y  UAB=0.0 

>  37.4  END 

>  37.5  PARAMETER  X=274,0y  IC1=6.0»  UAB=0.70 

>  37 . 7  END 

>  38  STOP 

>  39  ENDJOB 

♦END  OF  FILE 

# 


' 


